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ABSTRACT 


This  dissertation  describes  r  spec i al -purpose  siqnal 
processor  for  performing  two-d  imen  s  i  one  ]  convolution  with  a 
minimum  amount  of  hardware  using  the  concepts  of  sinqular 
value  decomposition  (BVD)  and  small  generating  kernel  (CCK) 
convolution.  The  SVD  of  an  impulse  response  of  a 
two-dimensional  finite  impulse  response  (FTPl  filler  is 
employed  to  decompose  a  filter  into  a  sum  of 
two-dimensional  separable  linear  operators.  The?r  linger 
operators  are  themselves,  decomoosed  into  a  sequence  of  small 
kernel  convolution  operators.  The  CVD  expansion  can  be 
truncated  to  a  relatively  few  terms  without  s iqn  i  f i can 1 1 y 
affecting  the  filter  output. 

A  statistical  analysis  of  finite  word-length  effects 
in  SVD/SGK  convolution  is  presented.  Two  important  issues, 
related  to  the  implementation  of  the  filters  in  cascade 
form,  scaling  and  section  ordering,  are  also  considered. 

Computer  simulation  of  image  convolution  indicates 
that  12  bits  are  required  for  the  BGK/SVD  accumulator 
memory  and  16  bits  are  required  for  quantization  of  filter 
coefficients  to  obtain  results  visually  indistinguishable 
from  full  precision  computation.  A  normalized  mean  souere 
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error  between  the  SVD/SGK  processed  output  end  the  direct 
processed  output  is  chosen  as  an  objective  criterion 
function.  It  is  shown  that  a  subjective  visual  improvement 
is  obtained  by  resetting  the  output  mean  to  be  eoual  to  the 
input  mean. 

The  tr ansformat ion  technique  developed  for  the 
one-dimensional  case  is  used  to  parametr  i  cal  1  y  modify  rhe 
cutoff  frequency  of  a  baseline  SVD/FC-K  convolution  filler. 
A  detailed  discussion  of  the  one-dimensional  case  is 
presented,  and  its  applicability  to  SVD/FGK  convolution 
filters  is  described. 
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CHAPTER  1 


INTRODUCTION 


During  the  last  decade,  the  field  of  digital  signal 
processing  has  been-  extremely  dynamic  and  active.  There 
have  been  many  applications  of  digital  signal  processinq 
techniques  in  digital  communication,  seismic  processing, 
radar  processing,  sonar  processing,  speech  processing,  and 
image  processing. 


One  of  the  important  areas  in  digital  signal 
processing  is  digital  filtering.  The  term  "digital 
filtering"  can  be  viewed  as  a  computational  process  or 
algorithm  by  which  a  sampled  signal  or  a  sequence  of 
numbers,  acting  as  an  input  signal,  is  transformed  into  a 
second  sequence  of  numbers  called  the  output.  There  are 
two  major  types  of  digital  filters:  infinite  impulse 
response  (HR)  filters  and  finite  impulse  response  (FIR) 
filters.  Digital  filtering  is  mainly  concerned  with  filter 
design  and  its  implementation. 


If  the 
the  past, 
inputs,  the 


present  output  of  a  system  is  calculated  from 
present,  and,  in  the  noncausal  case,  future 
system  is  called  nonrecursive.  If  the  present 


output  of  a  system  is  calculated  from  the  past  and  present 
inputs  and  outputs,  the  system  is  called  recursive.  In 
both  recursive  and  nonrecursive  systems,  the  relation 
between  an  input  sequence  x(n)  and  an  output  seauence  y(n) 
can  be  characterized  by  a  difference  ecuation  of  the  form 

M  N 

y(n)  =  ^  ^akf  (n-k) ~]bky  (n-k)  (1-1) 
k=0  k=0 

Conceptually,  M  and  N  can  be  finite  or  infinite.  A  system 
in  which  b^  =  0,  for  k  =  ],..., N,  is  nonrecursive,  and  can 
be  implemented  by  an  FIP  filter.  The  system  in  which  F>] 
and  bk  is  not  zero  is  recursive,  and  it  can  be  implemented 
by  an  HR  filter.  Choosing  between  an  FIR  filter  and  an 
IIR  filter  depends  upon  the  relative  advantages  and 
disadvantages  the  filter  offers  for  a  specific  problem. 

Signal  processing  is,  of  course,  not  limited  to 
one-dimension.  Nany  signals  are  inherently 
two-dimensional;  thus,  two-dimensional  signal  processing 
techniques  are  required.  Image  data  is  a  typical 
two-dimensional  signal.  Digital  filtering  with  FIR  filters 
has  many  applications  in  image  processing.  For  instance, 
image  restoration  to  remove  blur  and  to  suppress  noise 
generally  requires  digital  filtering.  In  most  cases, 
digital  filtering  requires  implementat ion  of  a 
two-dimensional  convolution. 
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The  term  "implementation"  means  that  the  algorithm  is 
either  written  in  a  computer  language  for  a  general-purpose 
computer  or  is  realized  with  special-purpose  hardware.  In 
general,  the  implementation  of  two-dimensional  convolution 
in  image  processing  has  been  confined  primarily  to  computer 
programs  with  a  general-purpose  computer,  where  virtually 
unlimited  memory,  processing  capability,  and  time,  are 
readily  available.  .But.  the  required  processing  time  is 
often  quite  enormous  because  of  the  huge  amounts  of  data  to 
be  processed  and  the  restricted  input-output  transfer  rate 
between  the  computer  and  display.  An  image  size  of  542x512 
pixels  is  common  in  image  processing.  An  alternative  to 
the  use  of  a  general-purpose  computer  is  to  utilize 
Integrated  Circuit  (TC)  technology.  Recent  advances  in  TC 
technology  now  make  the  realization  of  a  real-time  signal 
processor  capable  of  performing  two-dimensional  convolution 
practical.  High  speed  digital  multipliers,  memory  and 
display  circuitry  are  now  commercially  available.  As  a 
result,  significantly  more  sophisticated  algorithms  can  now 
be  chosen  for  problem  solving.  The  trend  is  to  develop 
special-purpose  signal  processors  to  take  advantage  of 
recent  developments  in  digital  circuits  [1-1  to  1-31.  In 
the  design  of  such  a  special-purpose  signal  processor, 
speed,  complexity,  power  consumption,  computing  capability, 
and  cost,  are  all  factors  to  be  considered. 
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Recently,  a  technique  called  small  generating  kernel 
{ SGK )  convolution  has  been  proposed  [  ] —  4 ]  .  SGK  convolution 
is  a  processing  technique  for  performing  convolution  on  a 
two-dimensional  data  array  by  sequentially  convolving  the 
array  with  a  small  size  convolution  kernel,  say  3x3.  This 
idea  was  first  suggested  by  Mersereau  et  al .  [1-5]  and 
generalized  later  by  Abramatic  and  Faugeras  [1-4].  Fincc-  a 
large  size  kernel  convolution  is  performed  by  a  sequential 
small  size  kernel  convolution,  and  the  implementation  is 
highly  modular,  the  SGK  approach  makes  the  hardware 
implementation  quite  appealing  if  a  proper  design  procedure 
to  specify  the  small  size  kernel  operators  is  found. 


In  the  one-dimensional  case,  any  impulse  response  can 
be  decomposed  into  small  size  convolution  operators, 
typically  3x1.  This  property  can  be  seen  in  the  cascade 
form  for  FIR  filters.  But,  theoretically,  exact 
decomposition  of  a  large  size  convolution  operator  into 
small  size  convolution  operators  is  impossible  in  the 
two-dimensional  case.  This  is  the  reason  why  the  design 
procedure  for  SGK  convolution  leads  to  a  complicated  and 
time-consuming  optimization  problem.  The  inherent 
difficulty  in  finding  small  size  convolution  operators 
motivates  the  development  of  a  new  algorithm  for  the 
two-dimensional  convolution.  The  proposed  SVD/SGK 
convolution  method  also  makes  use  of  SGK  convolution, 

however,  the  size  of  small  size  convolution  operators  is 
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3x1,  rather  than  3x3. 


This  dissertation  describes  a  special-purpose  signal 
processor  with  a  minimum  amount  of  hardware  for  performing 
two-dimensional  convolution  using  the  concepts  of  singular 
value  decomposition  (SVD)  and  SGK  convolution.  To  extend 
the  usefulness  of  SGK  convolution,  two-dimensional  FIR 
filters  of  size  N-^xN2  are  decomposed  into  a  sum  of 
two-dimensional  separable  filters  by  means  of  the  SVD  of 
their  impulse  response  matrix  H.  The  SVD  expansion  can  be 
truncated  to  K  terms  (K<  R,  where  R  is  a  rank  of  H)  , 
without  significantly  affecting  the  output  of  the  filter. 
Whenever  the  two-dimensional  FIR  filter  is  separable,  the 
convolution  can  be  performed  by  one-dimensional  processing. 
This  is  a  reason  why  the  SVD  expansion  can  be  very  useful 
for  implementing  two-dimensional  nonseparable  filters.  It 
was  noted  that  each  one-dimensional  FIR  filter  can  be 
realized  as  a  cascade  of  second-order  SGK  filters.  Thus, 
it  is  possible  to  implement  a  two-dimensional  convolution 
by  a  sequential  convolution  with  one-dimensional 
second-order  SGK  filters.  As  an  example,  one  can  think  of 
using  such  a  convolution  technique  for  convolving  images  at 
real-time  rates  on  an  image  display  system. 

When  a  digital  signal  processing  algorithm  is 
implemented  with  a  special-purpose  signal  processor, 
account  must  be  taken  of  the  errors  caused  by  finite 
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word-length  in  representing  filter  coefficients  and  signal 
values.  Implementation  with  finite  word-length  can  be 
modeled  by  injecting  white  noise  into  signals  whenever  a 
rounding  operation  is  performed.  The  goal  of  this  error 
analysis  is  to  minimize  the  required  word-length  subject  to 
some  reasonable  error  tolerance.  The  problem  is  to 
determine  the  best  ordering  and  scaling  procedure  in  order 
to  minimize  the  required  word-length.  To  solve  these  two 
problems,  we  show  that  how  the  theory,  for  the 
one-dimensional  case,  can  be  modified  to  the 
two-dimensional  case. 

The  second  issue  investigated  in  this  dissertation  is 

parametric  design.  The  concept  of  parametric  design  is  to 

generate  a  class  of  two-dimensional  FIR  filters  with  a 

variable  cutoff  frequency  from  previously  designed  baseline 

SVD/FGK  convolution  filters.  Tn  the  case  of 

one-dimensional  FIR  filters,  Oppenheim  et  al .  11—61  have 

proposed  a  transformat  ion  for  designing  a  variable  cutoff 

digital  filter.  But,  very  little  work  has  been  reported  in 

the  two-dimensional  case.  It  is  shown  that  the  cutoff 

frequency  of  a  SVD/FGK  convolution  filter  can  be  varied  by 

the  use  of  a  one-dimensional  transformation.  It  is 

believed  that  such  a  variable  cutoff  SVD/FGK  convolution 

filter  has  numerous  applications  in  image  processing. 

Adaptive  filtering  will  be  very  useful  in  image 

restoration.  For  example,  the  cutoff  frequency  of  a  Wiener 
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filter  could  be  changed,  and  an  observer  could  effectively 
examine  the  processed  image  in  real-time. 

This  dissertation  consists  of  seven  chapters.  A 
review  of  SCK  and  SVD/SCK  convolution  is  presented  in 
Chapter  2.  Chapter  3  discusses  the  effect  of  using 
fixed-point  arithmetic.  This  chapter  includes  a  derivation 
of  the  noise  formula  to  predict  total  roundoff  noise. 
Scaling  and  section  ordering  for  SVD/SGK  convolution  are 
described  in  Chapter  4.  In  Chapter  5,  a  series  of 
experimental  results  based  on  computer  simulation  is 
presented.  Among  these  results  is  the  confirmation  of  the 
derived  noise  formula  using  a  random  number  array  as  an 
input.  A  simple  technique  to  reduce  the  normalized  mean 
square  error  (NMSE)  between  the  SVD/SGK  processed  output 
and  the  direct,  processed  output  is  also  given.  The 
effectiveness  of  this  technique  is  demonstrated  visually. 
Chapter  6  deals  with  the  parametric  design  of  SVD/SGK 
convolution  filters.  A  detailed  discussion  of  the 
one-dimensional  case  is  presented,  and  its  applicability  to 
SVD/SGK  convolution  filters  is  described.  Several  design 
examples  for  two-dimensional,  as  well  as  one-dimensional 
cases,  are  shown  in  this  chapter.  rinally.  Chapter  7 
discusses  the  conclusion  and  possible  extension  of  this 
work. 
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CHAPTER  2 


SEQUENTIAL  CONVOLUTION  TECHNIQUES 

2 . 1  Introduction 

Two-dimensional  convolution  has  found  numerous 
applications  within  the  field  of  two-dimensional  signal 
processing  [2-1, 2-2].  For  example,  image  enhancement, 
image  restoration,  and  digital  filtering  generally  require 
two-dimensional  convolution.  Referring  to  Fig.  2-1,  an 
output  array  G(j,k)  is  obtained  by  convolving  the  input 
array  F(j,k)  with  the  impulse  response  of  the  system 
H(j,k).  The  two-dimensional  direct  convolution  algorithm 
can  be  expressed  by  the  double  sum 

j  k 

G ( j , k ) =F ( j , k ) »»H ( j , k )  =y  ^  ^F  (m,n)  H  ( j-m+1 ,  k-n+1)  (2-1) 

m=l  n=l 

where  G(j,k)  is  the  M^xf^  output  array,  F(j,k)  is  the  N^xN2 
input  array,  and  H(j,k)  is  the  L^xL2  convolution  kernel 
array,  called  an  impulse  response.  The  input  and  output 
dimensions  are  related  by  =  N^L^-l  for  i  =  1 , 2 .  In  Fq. 
(2-1),  the  symbol  ®@  denotes  a  two-dimensional 

convolution.  The  symbol  ®  will  be  used  to  represent  a 
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SjBF  nwn»i'W» '>»»»«*- 


one-dimensional  convolution  throughout  this  dissertation. 


In 

the 

direct  convolution 

algorithm, 

the.  output  , 

G(j ,k)  , 

is 

the  weighted  sum  of 

all 

vsl ues 

of  the  input 

array. 

The 

drawback  of  using 

the 

d  i  reef 

conv o 1 ut i on 

algorithm  of  Eci.  (2-1)  is  that  it  reauires  many  arithmetic 
operations.  The  number  of  additions  and  multiplications 
for  direct  convolution  is  ^3^2  L'l  L'2  • 

In  1965,  Turkey  and  Cooley  [2-31  opened  a  new  era  in 
digital  signal  processing.  They  discovered  a  fast  Fourier 
transform  (FFT)  algorithm,  which  is  an  efficient  method  for 
computing  a  discrete  version  of  the  Fourier  transform 
(DFT) .  The  two-dimensional  DF1  pair  of  a  finite  array 

X ( j , k )  for  j , k  =  0,1 . N-l  can  be  written  in  the  form 

N-l  N-l 

X (u, v) =  -  SEE  X(j,k)exp|-  ?~(uj+vk)j- 
N  j=0  k=0 

(2-2) 

N-l  N-l 

X  ( j  ,  k )  =  EE  X  (u ,  v)  exp|“i  (u  j+vk  )| 
u=0  v=0 

where  i  =  VT.  u,v  are  spatial  frequency  variables,  and 
X  (u ,v)  denotes  the  Fourier  transform.  Both  %(u,v)  and 
X(j,k)  are,  in  general,  complex  series.  Consider  the 
following  relation  in  the  frequency  domain 


ir(u,v)  =  3(u,v)V(u,v) 


(2-3) 
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where  ^(u,v),  <i(u,v)  and  V(u,v)  are  discrete 


Fourier 


transform  of  the  array  G(j,k),  F(j,k),  and  H(j,k), 

respectively.  By  the  definition  of  the  DFT,  C(j,k)  can  be 
expressed  as 

Nl-1  N2“1 

G(^'k)=5Z  2  [/*(u,v)*  (u,v)jexp|2-i  (p-  +  (2_4) 

u=0  v=0  1  2 

Thus,  computation  of  the  discrete  convolution  of  two  arrays 

can  be  obtained  indirectly  using  the  DFT.  Considerable 

computational  efficiency  can  be  gained  by  the  FFT 

convolution  technique.  In  general,  computation  requires 
2  2 

N  +4N  log2  N  operations  when  N  [2-4]. 

Fourier  domain  processing  is  more  computat ional 1 y 
efficient  than  the  direct  convolution  of  Eq.  (2-1)  if  the 
size  of  the  impulse  response  is  sufficiently  large.  The 
cross  over  point  for  the  two  implementations  occurs  for  a 
10x10  impulse  response  with  large  input  arrays  [2-5], 
Because,  in  many  practical  applications,  the  size  of  an 
impulse  response  is  larger  than  10x10,  then  Fourier  domain 
processing  is  an  efficient  computation  technique. 
Furthermore,  the  efficiency  of  Fourier  domain  processing 
can  be  increased  by  overlap-add  or  overlap-save  techniques 
[2-6]  . 


Several  other  techniques,  for  example,  number 
theoretic  transforms,  have  been  reported  concerning 
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convolution  computation  [2-7, 2-8],  So  far,  the  technioues 
we  have  discussed  can  be  implemented  by  proarems  for  a 
general-purpose  computer  or  special-purpose  hardware. 
Recently,  due  to  the  dramatic  development  in  Large  Scale- 
Integrated  (LSI)  circuit  technology,  real-time  low  cost 
hardware  implementation  of  a  two-d i mens ional  convolution  is 
of  great  interest.  Low  cost  hardware  implementation  is 
possible  if  the  size  of  the  convolution  kernels  is  kept 
small  because  the  cost  of  hardware  is  proportional  to  the 
size  of  the  convolution  kernel.  The  techniauo,  commonly 
referred  to  as  PGK  convolution,  makes  this  task  possible. 
A  review  of  these  methods  is  given  in  Section  2-2.  The 
basic  concepts  of  the  SVD  technioue  dealing  with 
nonseparable  impulse  response  and  application  to  sequential 
convolution  is  discussed  in  Section  2-3.  A  new  convolution 
technique  is  proposed  in  Section  2-4.  Its  application  to 
an  image  processing  display  system  is  described  in  Section 
2-5. 

2 . 2  Review  of  Small  Generating  Kernel  Convol ut ion 

SGK  convolution  is  a  processing  technioue  for 
performing  convolution  on  two-dimensional  data  arrays  by 
sequentially  convolving  the  arrays  with  small  size 
convolution  kernels.  The  output  of  the  SGK  convolution 
operation  closely  approximates  the  output  obtained  by 
convolution  with  a  large  kernel  prototype  filter.  The 
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motivation  behind  FCK  convolution  is  that  it  can  be  used  to 
approximate  any  impulse  response  of  an  FT  R  filter,  and  that 
its  structure  permits  implementation  of  the  convolution  by 
sequential  convolution  with  small  size  kernels. 

McClellan  12-S]  was  the  first  to  propose  a  techniaue 
for  designing  such  a  class  of  filters  by  transforming 
one-dimensional  linear  phase  filters*  into  two-dimensional 
linear  phase  filters.  Ey  assuming  that  the  prototype 
filter  is  a  linear  phase  filter,  his  algorithm  transforms  a 
one-dimensional  filter  Mu)  into  a  two-dimensional  filter 
%'  (u,v)  by  means  of  transformation  given  by 

cosw  =  Acosu+Bcosv+Ccosu -cosv+D  (2-5) 

The  McClellan  transformation  is  an  extremely  useful  tool, 
requiring  only  moderate  computation,  for  designing  many 
common  types  of  two-dimensional  FIR  filters.  FIR  filters 
up  to  order  100  can  be  designed  using  this  method. 

Mersereau  et  al  .  [2-10]  generalized  the  McClellan 
transformation  for  two-dimensional  FIR  filters  and  showed 
an  efficient  way  to  implement  the  filters  designed  by  this 
method.  The  significance  of  their  implementation  of  the 
designed  filter  is  that  a  large  two-dimensional  convolution 


*A  linear  phase  filter  implies  symmetry  of  the  filter. 
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can  be  replaced  by  a  sequential  convolution  with  small  size 
kernel  operators.  A  description  of  the  algorithm  follows. 


The  frequency  response  of  a  one-dimensional  linear 
phase  filter  of  odd  length  L  is 


L-l 

2 

Mu)  =  ^>  'h  (n)  [cosuln 
n=0 


(2-6) 


where  h(n)  represents  the  filter  impulse  response.  Because 
the  frequency  response  of  a  cascade  form  is  the  product  of 
the  frequency  response  of  individual  stages,  the  term 
[cos  u]n  of  Eq.  (2-6)  can  be  considered  as  a  total 
frequency  response  obtained  by  cascading  n  identical 

filters  each  with  a  frequency  response  cos  u.  It  is 

beneficial  to  rewrite  Eq.  (2-6)  in  terms  of  the  z-transform 
to  obtain 

L-l 

L-l  2 

H  (z)  =  (n)  z  n=h(0)+^~\(n)  [p1  (z)  ]n  (7-1  s) 

n=0  n=0 


where 


px  (z: 


z+z 


-1 


( 2  —  7  b ) 


Figures  2-2  to  2-4  show  three  basic  implementation 
structures  proposed  by  Mersereau  et  al  .  [2-111. 
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Refecting  to  Fig.  2-2,  implementation  of  a 

L-l 

two-dimensional  filter  consists  of  a  Q  (=  ~2~  )  stage 
identical  sequential  convolution.  Note  that  p^(z)  is 
replaced  by  a  two-dimensional  filter  Ff(z  ,z2),  which  is 
obtained  by  the  NcClellan  transformation.  The  q-th  staqe 
output  Cq(Zj,Z2)  is  obtained  from  the  cumulative  sum  of  the 
q-th  stage  as 

VW  =  Oq_1(z1/z2)+h(q)Aq(z1,Z2)  (2-8?) 

wh  e  r  e 

Aq(2l'Z2)  =  Aq-l(zl'Z2)Hf (Z1'Z2)  f  2  —  8  b ) 

The  term  Oq(z^,Z2>  corresponds  to  the  output  array 
G(z1,z2)/  or  equivalently 

Q 

°q(z1,z2>  =  ^h(U(Hf(z1,z2)]JF(z1,z2)  (2-9) 

t=l 


The  convolution  indicated  in  Eqs.  (2-8)  and  (2-9)  could  be 
implemented  directly  from  the  direct  convolution  algorithm 
of  Ea.  (2-1).  The  other  structures,  shown  in  Figures  2-3 
and  2-4,  also  can  be  implemented  in  a  similar  manner 
[2-11].  Mersereau  et  al .  also  pointed  out  that  the 
computational  efficiency,  i.e.,  number  of  multiplication 
and  addition  operations  required  for  implementation,  is 
greater  than  the  number  for 


either 


the  direct 
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F(z,,z2) 


implementation  ot  an  FFT  impl ementat ion  for  filters  of 
order  up  to  50x50. 

Another  structure  of  interest,  shown  in  Fig.  2-5,  was 

proposed  by  McClellan  and  Chan  [2-12].  They  noted  that 
n  -n 

_  is  the  n-th  degree  Chebyshev  polynomial  in  the 

variable  p^(z)  of  Eq.  (2-6).  Unfortunately,  an  arbitrary 
two-dimensional  impulse  response  cannot  be  implemented  in 
this  way  because  it  is  not  always  symmetrical.  The 
implementations  discussed  so  far  are  applicable  only  to 
McClellan  transformed  filters.  The  elementary  filters  of 
Figures  2-2  to  2-4  do  not  necessarily  has  the  same 
frequency  response.  The  limitation  of  the  previous 
implementation  has  motivated  a  search  for  more  general 
design  techniques  for  a  class  of  two-dimensional  FIR 
filters  that  can  be  easily  implemented  by  sequential 
convolution  with  small  size  kernels,  say  3x3. 

Abramatic  and  Faugeras  f  2—1 3  to  2-15]  presented  a 
synthesis  procedure,  described  in  Fig.  2-6,  for  designing 
such  a  class  of  filters.  In  comparison  with  Fig.  2-2,  the 
elementary  second-order  filters  have  different  transfer 
functions.  The  sequential  filter  proposed  by  Mersereau  et 
al  .  is  a  special  case  of  this  class  of  filters.  The  design 
procedure  approximates  the  prototype  filter  by  means  of 
minimizing  the  mean  square  error  [2-13]  or  Chebyshev  error 
[2-14]  between  the  approximate  and  prototype  filters. 
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H  ( 1 , 2  J 

H(2,l) 


H(1,L2) 


H<L1,L2! 


(2- 30) 


Suppose  the  impulse  response  is  spatial  i  • 

of  separable  f  ^  lpvariant  and  jc 

Pcrrb]f  form  such  that 
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(2-1  1  ) 


where  c  and  £  are  column  vectors  representing  column  and 
row  one-dimensional  impulse  responses,  respectively.  We 
have  used  the  superscript  T  to  denote  transposition.  Then, 
two-dimensional  convolution  may  be  performed  by  sequential 
row  and  column  one-dimensional  convolutions.  As  a  result, 
one  can  obtain  a  substantial  decease  in  the  number  of 
multiplication  and  addition  operations  if  the  input  array 
size  becomes  large.  If  the  input  array  size  is  NxN,  the 

separable  convolution  operators  of  Eq.  (2-11)  reauires 

2  .  .  2 
N  (L^+L2)  multiplications  compared  with  N  L2 

multiplications  required  in  the  nonseparable  case  (fewer 

are  required  if  the  impulse  response  matrix  possesses 

symmetry)  .  Unfortunately,  we  cannot  assume  that  the 

prototype  impulse  response  matrix  H  is  always  separable. 

One  way  of  dealing  with  the  nonseparability  problem  is  to 

use  the  SVD  technique.  In  the  SVD  matrix  expansion,  any 

arbitrary  L^xL^  matrix  of  rank  R  can  be  decomposed  into  the 

sum  of  a  weighted  set  of  unit  rank  L-^xI^  matrices.  The 

significance  of  the  SVD  expansion  is  demonstrated  by  noting 

that  the  nonseparable  matrix  H  is  the  sum  of  individual 

separable  matrices  of  unit  rank  [2-17]. 

According  to  the  SVD  expansion,  there  exist  an  L-^xL^ 
unitary  matrix  U  and  an  I^xL^  unitary  matrix  V  for  which 

T 

U  H  V  =  A 2 


(2-12) 
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where 


3, 


■x'£d)  0 

0  " 

A%(2) 

0  'A^(R) 

— 

0 

0 

}, 

is  an  L^xL^  matrix  with  a  general  diagonal  entry  -  (j)  for 
j=l,2,...F  called  a  singular  value  of  H.  The  singular 

values  can  be  obtained  by  square  rooting  the  eigenvalues  of 

T  T  T 

HH  or  H  H.  The  columns  of  U  are  the  eigenvectors  of  HH 

T 

and  the  columns  of  V  are  the  eigenvectors  of  H  F.  Since 
T  T 

HH  and  H  H  are  symmetrical  and  square,  the  eigenvalues 
X^(j)  are  real,  and  the  eigenvectors  set  {u^}  •  *  v  j  )  for 

j  =  1,2, ...R  are  orthogonal. 


Since  matrices  U  and  V  are  unitary  matrices, 
Eq.  (2-12)  is  equivalent  to  Eq.  (2-14).  Hence,  H  can  be 
decomposed  as 


H  =  (Uj  u2  . 


(2-14) 


Equation  (2-14)  can  be  reformulated  into  series  form  as 
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(2-15) 


If  we  let 


S)  * 


(2-1 6a) 


r  .  =  v  . 
.  “3  ~3 


2-1  6b) 


where  c  .  and  r  .  are  one-dimensional  column  and  row 
-  1  —  3 

convolution  operators,  respectively,  then 


( 2  —  1 7  a ) 


where 


H  .  =  c  .  •  r  . 
-3  -3  -3 


( 2-1 7b) 


It  should  be  observed  that  the  vector  outer  product  u.»v. 

_3  --j 

of  the  eigenvectors  forms  a  set  of  separable  unit  rank 

matrices  each  of  which  is  weighted  by  a  corresponding 

singular  value  of  H,  as  shown  in  Fig.  2-8.  If  matrix  H  is 

separable,  then  we  have  only  one  SVD  expansion  term.  If 

matrix  H  is  not  separable,  theoretically,  the  exact 

representation  of  H  needs  R  terms.  Hence,  the  number  of 

multiplication  operations  for  direct  convolution  requires 

RN2(L1+L2)  multiplication  operations,  as  shown  in  Fig.  2-9. 
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Figure  2-8.  SVD  expansion  into  unit  rank  matrices 


are 


If  we  assume  that  the  singular  values  )S  (j) 
listed  in  order  of  decreasing  magnitude,  then  the  SVD 
expansion  of  Eq.  (2-17a)  can  be  always  rewritten  as 


(2-1  8b) 

(2-1 8c) 


where  K  is  the  number  of  retained  term  for  HSVD  and  the 
"Hat"  symbol  (  ^  )  represents  the  approximation  of  H.  The 
matrix  ER  denotes  the  truncation  error  as  a  result  of 
retaining  the  first  K  terms.  Obviously,  E  =0  for  K  =  R. 

-  K 

It  can  be  shown  that  the  case  for  K  =  1  corresponds  to  the 

minimum  mean  square  error  (NMSE)  separable  approximation  of 

H  [2-18].  If  the  elements  of  H  for  K=1  are  not  close 

“SVD 

in  comparison  to  the  elements  of  H,  we  may  add  the  next 
largest  singular  value  term  for  a  closer  approximation  of 


H. 


In  general,  we  will  be  satisfied  with  a  multi-stage 
expansion  that  will  closely  approximate  H.  One  of  the  most 
commonly  used  numerical  error  measurements  is  the 


normalized  mean  square  error  (NMSE).  Let  us  define  the  SVD 
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approximation  error 

eK  for  the  measurement  of 

the 

degree 

of  approximation  by 

retaining  only  the  first  K 

terms 

in  the 

expansion  as 

EE'v^’i2 

ESH<i'3'i2 


(2-19) 
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In  most  cases,  an  imaging  system  can  be  modeled  by  a 
superposition  integral  relating  the  input  and  output 
continuous  fields  of  a  linear  system  [2-4].  In  order  to 
reduce  the  problem  to  a  discrete  model,  the  point  spread 
function  (PSF)  of  the  imaging  system,  as  well  as  the  input 
and  output  images,  should  be  discretized.  The  matrix  H 
resulting  from  the  PSF  samples  is  nearly  singular  or 
ill-conditioned  since  the  rows  of  the  matrix  H  are 
approximately  a  linear  combination  of  one  another 
[2-19,2-20] . 
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Ill-conditioning  of  a  matrix  can  be  described  by  its 
condition  number  [2-21].  The  condition  number  of  given 
matrix  A  is  defined  in  terms  of  the  ratio 


C  [A] 


max 


mm 


(2-20) 


"2  'V 

of  the  largest  \  to  smallest  \ '  singular 

max  mm 

value  of  A.  The  condition  number  is  a  useful  tool  for 
explaining  the  effect  of  perturbation  caused  by  additive 
noise  on  the  accuracy  of  computation  involved  [2-22].  The 

i, 

condition  number  approaches  infinity  as  A2-n  goes  to 

zero.  In  this  case,  the  matrix  is  called  ill-conditioned 
and  will  have  a  large  condition  number.  In  an  idea] 
imaging  system,  characterized  by  a  delta  function  point 
spread  function,  the  condition  number  is  unity  since  all 
singular  values  have  the  same  magnitude.  Sometimes  it  is 
convenient  to  demonstrate  matrix  conditioning  by  showing 
singular  value  magnitude  plots.  Referring  to  Fig.  2-10,  a 
well-conditioned  matrix  requires  more  terms  in  a  SVD 
expansion  than  an  ill-conditioned  matrix.  But,  it  is  noted 
here  that  ill-conditioned  and  nearly  singular  problems  are 
very  common  in  imaging  systems  [2-4].  Therefore,  we  do  not 
need  to  retain  all  terms  in  the  SVD  expansion,  but  only  a 
few  terms  because  of  ill-conditioning  of  the  PSF  matrix 
itself.  The  usefulness  of  the  SVD  expansion  can  be 
demonstrated  by  noting  that  we  can  trade  off  between  the 
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a)  ill -conditioning 


b)  good -conditioning 


c)  ideal- conditioning 


Figure  2-10.  Singular  value  plot 


amount  of  NMSE  and  the  computational  efficiency  by  choosing 


the  number  of 
only  K  terms 
multiplication 
still  holds  as 


terms  in  the  SVD  expansion.  By 
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2 
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2 • ^  The  SVD/SGK  Cascade  Convol ut ion  Techn ique 


In  the  previous  section,  approximion  of  a  nonsep-rable 
impulse  response  matrix  H  in  terms  of  the  sum  of  individual 
separable  matrices  of  unit  rank  was  discussed.  To 
implement  the  SVD  convolution,  each  separable  convolution 
operator  is  implemented  in  parallel,  and  summed  together, 
as  shown  in  Fig.  2-9.  In  this  section,  SVD  and  SGK 
techniques  are  combined  to  obtain  a  more  versatile 
two-dimensional  convolution  technique  requiring  a  simpler 
implementat  ion . 


Since  each  SVD  expanded  separable  matrix  of  unit  rank 

is  an  outer  product  of  the  one-dimensional  column  and  the 

row  operator  Cj  and  _r  j,  here  each  Cj  and  £j  is  to  be 

considered  as  a  one-dimensional  FIR  filter.  There  are  a 

variety  of  alternative  forms  -in  the  FIR  filter  realization. 

Realization  of  FIR  filters  generally  takes  the  form  of  a 

nonrecursive  computation  algorithm.  One  way  of  realizing 

FIR  filters  for  hardware  simplicity  is  to  use  a  cascade 

form.  In  the  cascade  form,  the  z-transform  of  the  impulse 

response  with  the  length  of  L  can  be  expressed  as  a  product 
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of  second-order  SGK  filters  as 


Q  Q 

H(z)=  nHk(z)=  n  I60(k+81/kz'1+62ikz'2] 

k=l  k=l 


(2-2]  ) 


where  the  8.  .  are  real  numbers  and  Q,  the  number  of 

j  i K 

convolution  stages,  is 


Q  = 


L-l 

2 

L 

2 


,  L  odd 
,  L  even 


(2-22) 


When  L  is  even,  one  of  the  kernel  terms  B2  jcwill  be  zero. 
Here  we  shall  be  concerned  only  with  the  case  of  odd  length 
impulse  response.  The  kernel  of  each  second-order  SGK 
filter  can  be  easily  obtained  by  solving  the  zeros  of  the 
polynomial  H(z)  because  H(z)  is  a  polynomial  in  z  *  of 
degree  L-l . 

A  new  approach  for  two-dimensional  SVD/FGK 

convolution,  shown  in  Fig.  2-11,  is  to  realize  each 

one-dimensional  convolution  operator  c  and  r  for 

-  £  -  £ 

£  =  1,2,... , K  as  a  sequence  of  second-order  SGK  filters. 
Referring  to  Fig.  2-11,  the  z-transform  of  the  SVD/SGK 
convolution  filter  is 

K 

H(Zi'Z2)  -  (z^R^  ( z2 5  (2-233) 

£=1 

or 
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( 2  —  2  3  b ) 


'•■h  ere 


Q 

C f.  ( z | )  —  n  c,  ■  (Z  )  (2  —  243) 

J-  i=1  x 


Rp(z2) 


Q 

n  r,  ,  (z 2> 

i=l  -'J  c 


( 2-2  4b) 


The  terms  and  P£(z2)  for  £=1,2,. ..K  denote  the 

2-transform  of  each  column  and  row  one  dimensional 
convolution  operator,  as  defined  in  Eq.  (2-16),  and  each 
Cj  f(  z  i>  ,  R{  ^(z2)  for  i,j  =  1 ,  .  .  .  C  is  the  z-transform  of  the 
second-order  SGK  filter. 


One  of  the  most  important  reasons  for  using  FIR 
filters  is  that  they  can  be  designed  to  possess  linear 
phase,  a  feature  that  is  very  useful  in  speech  processing 
and  data  transmission.  It  is  easy  to  see  where  the  zeros 
of  such  linear  phase  FIR  filters  can  lie  by  examining  their 
z-transforms  because  a  linear  phase  filter  is  symmetrical. 
In  the  general  case,  the  filter  system  function  is 

L-l 

H<2)  =Eh(n,Z'n  (2-25) 

n=0 


Linear  phase  FIR  filters  have  a  symmetry  property  such  that 
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h(n)  =  h(L-l+n) 


(2-26) 


Therefore,  by  using  Eq.  (2-26),  Eq.  (2-25)  can  be  rewritten 


(L-l) 
2  ", 


H  ( z  )  =  z 


L-l  _  (L-l) 

h  ( 0 )  [z  2  +  z  2  1 

L-3  _  (L-3) 

+h  ( 1 )  }  [jz  2  +  z  2  J 


|  h(0)  [: 


(2-27 ) 


If  z  is  replaced  by  z  ,  then  we  obtain 


H  (z)  =  z 


(L-l) 

2  ( 


(L-l) 

2 


(L-l) 

2 


|  h  ( 0 )  C  z  2  +  z  2  ] 

1  -  (L-3)  (L-3) 

+h(l)  [z  2  +  z  2  ] 


(2-28) 


+  . 


By  comparing  Eq.  (2-21)  with  Eq.  (2-28),  one  obtains 


-(L-l) 


H  (z)  =  z 


H (z_1 ) 


(2-29) 


Equation  (2-27)  shows  that  the  zeros  of  H(z)  are  identical 


to  the  zeros  of  H(z" 


In  other  words,  if  H(z)  has  a 


_2  .  .2 


complex  zero  a+ib,  with  a  +b  /  1 ,  then  B(z)  must  have  a 

minor  image  zero  — .  Since  the  impulse  response  of 

a+ib 

the  filter  is  a  real  number,  every  complex  zero  of  H(z)  has 
its  complex  conjugate  as  another  zero. 


The  discussion  above  leads  immediately  to  the  possible 

2  2 

form  of  H,(z).  For  every  complex  zero  of  H(z),  a  +b  /  1, 
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the  second-order  SGK  filter  will  be  of  the  form 

Hk(z)  =  [z"1- (ak+ibk) ]  [z*1- (ak-ibk) ]  (2-30) 

If  the  zero,  a  ,  is  not  complex,  then  the  form  of  SGK 
filter  is 

H.  (z)  =  [z-1--a ]  [ z ~ 1  —  \]  (2-31) 

k  a 

If  the  zeros  are  either  -1  or  1,  then  the  zero  is  its  own 
complex  conjugate  as  well  as  s  mirror  image.  In  this  case, 
the  form  of  the  SGK  filter  is 

Hr(z)  =  (z_1t  1)  (2-32) 

From  the  discussion  above,  the  following  rule  of  zero 
grouping  can  be  stated : 

1)  Complex  zeros  are  grouped  together  in  conjugate 
pairs. 

2)  Heal  zeros,  that  are  reciprocal  of  each  other, 
are  paired  together. 

3)  Double  or  higher  multiplicity  zeros  are  paired 
together  in  pairs. 

The  rule  of  zero  grouping  guarantees  that  all  kernels 

are  real  numbers.  The  proposed  SVD/SGK  convolution  has 

both  advantages  and  d  isadvantages .  Since  two-dimensional 

large  kernel  convolution  is  replaced  by  a  cascade  of 

one-dimensional  SGK  filters,  the  processing  complexity  can 

39 


be  reduced.  Also,  from  a  theoretical  point  of  view,  there 
is  no  approximation  error  in  realizing  the  cascade  form 
because  all  kernels  can  be  found  exactly  by  simply  solving 
for  the  zeros  of  H(z).  Only  the  SVD  truncation  error 
defined  in  Section  2-3  will  be  introduced.  On  the  other 
hand,  computational  inefficiency  could  be  one  of  the 
disadvantages  of  replacing  two-dimensional  SGK  filters  by 
one-dimensional  SGK  filters'.  It  is  possible,  however,  to 
perform  two-dimensional  SGK  filter  convolution  instead  of 
one-dimensional  since  we  can  rewrite  Eq.  (2-23)  in  the 
alternative  form 


K 


where 


i<Vz2: 


’  Ct,i(zl'z2)Rt.i(2l'22l 


( 2-3  3b) 


As  a  matter  of  fact,  the  two-dimensional  SGK  filter  will 
increase  computational  speed  by  a  factor  of  two,  but  the 
hardware  is  more  costly  and  the  processing  more  complex. 
Implementation  of  a  two-dimensional  SGK  filter,  in  general, 
requires  nine  multipliers  and  adders. 

2 • ^  Image  Processing  Display  Impl ementat ion 

There  are  many  ways  to  implement  the  SVD/SGK 


convolution  method. 


The  goal  of  this  section  is  to 


describe  how  to  organize  the  implementation  and  apply 
SVD/SGK  convolution  to  an  image  processing  display  system. 
Let  us  denote  F(j,k)  as  a  filter  input  array  with  a  size  of 
NxN  and  the  array  G(j,k)  as  its  output.  We  also  assume 


that 

the  size 

of  the  prototype  impulse 

response 

is 

<2o+i : 

lx (2Q+1) . 

For  simplicity,  we  shall 

discuss 

only 

implementation  for  one  term  in  the  SVD  expansion  because 
the  SVD/SGK  convolution  consists  of  K  identical  paths.  The 
implementation  iterates  2Q  stages.  The  node  labeled 
Y i ( j  » k )  for  i  =  1 , 2 , . . . , ?Q  is  an  intermediate  array,  which 
will  be  used  in  the  next  convolution.  In  other  words,  at 
each  node,  the  array  ,k)  is  used  to  produce  array 

Yi(j,k).  Therefore,  Y2Q(j,k)  corresponds  to  the  final 
output  array  G(j,k).  Once  the  array  Yi(j,k)  has  been 

computed  from  Y  1(j,k),  Yi  ^(j,k)  is  no  longer  needed  and 

that  Y  ^  can  then  be  stored  in  place  of  Y^_^.  To  implement 
SVD/SGK  convolution,  it  is  then  necessary  to  have  at  least 
one  common  storage  for  the  intermediate  array  Y^(j,k)  for 
i  =  1,2,. ..,20.  Eut  the  storage  array  should  be 
initialized  by  the  input  array  F(j,k).  As  the  computations 
proceed  along  the  chain  of  SGK  filters,  each  Y^(j,k)  will 
be  larger  in  extent  than  its  predecessor  Y^  ^(j.k). 

Therefore,  the  required  storage  size  must  be  large  enough 
to  hold  (N+2Q) x (N+2Q)  pixels. 

Because  the  implementation  of  SVD/SGK  convolution  is 

highly  modular,  the  concept  of  SVD/SGK  convolution  is 
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ideally  suited  for  implementation  by  a  diqita]  imaqe 
display  system.  Two-dimensional  convolution  performed  by  a 
digital  computer  in  image  processing  is  often  quite  time 
consuming  because  of  the  serial  nature  of  the  computation 
and  the  slow  input-output  transfer  rate  between  the 
computer  and  display  [2-23].  But  solid  state  device 
technology  makes  it  possible  to  develop  memory  devices  that 
produce  pixels  at  a  serial  rate  of  about  10  million  per 
second.  Figure  2-12  is  a  basic  diagram  of  the  architecture 
for  SVD/SGK  convolution  [2-23].  Tn  the  operation  of  this 
hardware,  an  original  image  to  be  convolved  is  written  into 
an  accumulator  memory  with  a  size  of  (N+20) x (N+20) .  The 
accumulator  will  thus  appear  as  an  array  of  nonzero  values 
encircled  with  0  square  rings  of  zeros.  Then  the  input 
array  is  sequentially  convolved  with  a  3x1  impulse  response 
operator,  depending  on  the  row'  or  column  direction.  Three 
multiplication  and  three  addition  operations  are  performed 
for  each  pixel.  After  each,  convolution,  the  microprocessor 
will  update  the  kernels  of  the  3x1  convolution  operator. 
This  process  proceeds  for  20  stages,  equivelent  to  20  frame 
time.  Thus,  after  20  frame  times,  the  contents  of  the 
accumulator  memory  are  added  to  the  partial  sum  memory, 
which  is  initialized  by  zero,  and  return  to  the  original 
image.  This  processing  completes  the  first  term  of  the  BVD 
expansion.  The  partial  sum  memory  can  be  displayed,  if 
desired.  This  process  is  repeated  for  the  remaining  FVD 
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terms,  resulting  in  a  total  processing  time  of  2KQ  frame 
time  intervals.  For  conventional  30  f  r  ame/second 
operation,  the  SVD/SGK  convolution  operation  can  be 
completed  in  2KC/30  seconds,  far  less  than  the  20  to  30 
seconds  required  by  a  hardware  floating  point  processor. 

2 . 6  Conclusion 

In  this  chapter,  it  was  shown  that  the  SVD  expansion 
of  the  impulse  response  of  a  two-dimensional  FIR  filter  is 
a  very  useful  technique  for  a  two-dimensional  convolution. 
The  SGK  and  SVD/SGK  convolution  methods  are  attractive 
techniques  for  simplifying  the  computational  requirement  of 
two-dimensional  convolution.  The  SVD/SGK  convolution 
approach  is  attractive  for  two  reasons.  First,  large 
two-dimensional  convolution  is  replaced  by  seauential 
one-dimensional  convolution  with  small  size  convolution 
operators.  If  one  is  interested  in  implementing  SVD/SGK 
convolution  with  special-purpose  hardware,  that  approach 
reduces  both  the  cost  and  the  complexity  of  the  processing. 
Second,  the  design  for  the  SVD/SGK  convolution  filter  is 
simple  and  fast,  and  the  design  procedure  introduces  a  very 
small  approximation  error  caused  by  retaining  only  the 
first  few  terms  of  the  SVD  expansion.  But  the  design  for 
the  SGK  convolution  filter  generally  leads  to  complicated, 
time-consuming  nonlinear  optimization  programs.  To  utilize 
SVD/SGK  convolution  on  the  digital  image  display  system,  a 


bSSiC  dia5ra"  Of  architecture  to 
was  introduced. 


r  pvd/^GK  convolution 


CHAPTER  3 


FINITE  WORD-LENGTH  EFFECTS  IN  SVD/SGK  CONVOLUTION 

3  •  1  Introduction 

Until  now,  we  have  assumed  full  precision 
implementation  for  SVD/SGK  convolution.  We  will  now 
discuss  some  practical  problems  that  must  be  considered 
when  digital  signal  processing  algorithms  are  implemented 
with  programs  for  general-purpose  computers  or,  especially, 
with  special-purpose  hardware.  These  problems  are  caused 
by  the  use  of  finite  word-length  registers  to  represent 
signal  values,  coefficient  values,  and  arithmetic 
operations.  Because  of  finite  word-length,  a  quantized 
number  will  not  take  the  exact  value  assigned  by  the  design 
procedure . 

When  a  signal,  to  be  processed  digitally,  is  obtained 
by  sampling  a  band-limited  signal,  the  numbers  must  be 
represented  by  a  finite  word-length  register  in  the  digital 
machine.  This  conversion  process  between  analog  samples 
and  discrete  valued  samples  is  called  the  quantization 
process.  This  quantization  process  is  an  irreversible 
nonlinear  operation.  When  the  filter  coefficients  are 


quantized  for  digital  implementation,  the  resulting  filter 
must  be  checked  to  be  sure  that  it  is  close  enough  to  the 
desired  response.  In  addition,  finite  word-length 
operation  has  a  strong  effect  on  both  the  cost  and  speed  of 
the  system.  If  the  word-length  is  large,  then  the  cost  of 
hardware  will  be  expensive  and  the  processing  speed  low. 
Therefore,  reducing  the  w'ord-length  as  much  as  possible  is 
a  major  goal  . 

It  should  be  noted  here  that  effects  of  finite 
word-length  in  a  digital  filter  depend  on  many  issues  such 
as  whether  fixed-point  or  floating-point  arithmetic  is 
used,  whether  the  fixed-point  number  represents  a  fraction 
or  an  integer,  and  whether  quantization  is  performed  by 
rounding  or  truncating. 

In  a  digital  system,  numbers,  generally,  are 

represented  by  a  radix  of  2.  Thus,  numbers  are  represented 
by  strings  of  binary  digits,  either  zero  or  one.  If  a 
word-length  of  b  bits  is  chosen  to  represent  numbers,  2 
different  numbers  can  be  represented.  There  are  two  ways 
to  represent  binary  numbers,  depending  on  the  location  of 
binary  points.  In  fixed-point  arithmetic,  the  position  of 
a  binary  point  is  assumed  to  be  fixed.  The  bits  to  the 
right  of  a  binary  point  are  the  fractional  part  and  those 
to  the  left  of  the  binary  point  ere  the  integer  part.  But, 
with  no  loss  of  generality,  we  assume  throughout  this 
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dissertation  that  the  position  of  the  binary  point  is  just 
to  the  right  of  the  first  bit.  Thus,  the  range  of  numbers 
that  can  be  represented  with  b  bits  is  -1.0  to  1.0-2  ^ 

It  is  noted  that  the  signal  magnitude  can  be  scaled  to  any 
desired  range.  Certainly,  the  binary  point  could  be  moved 
further  to  the  right  to  allow  a  signal  with  magnitude 
greater  than  unity,  but  the  price  paid  is  greater 
complexity  in  hardware. 

There  are  three  formats  commonly  used  to  represent 
fixed-point  numbers,  depending  on  the  way  of  expressing 
negative  numbers:  sign  and  magnitude,  2's  complement,  l's 
complement.  The  sign  and  magnitude,  the  most  simple 
format,  represents  the  magnitude  by  a  binary  number;  the 
sign  is  represented  by  the  leading  digit.  It  is  useful  to 
assume  that  in  all  three  representations,  the  leading  bit 
is  zero  for  a  positive  number  and  one  for  a  negative 
number.  For  this  reason,  the  leading  bit  is  called  a  sign 
bit.  But  the  sign  and  magnitude  format  presents  an 
inherent  problem  in  performing  simple  arithmetic,  such  as 
addition.  Therefore,  the  sign  and  magnitude  format  is 
generally  avoided  in  a  digital  system. 

For  l's  complement  representation,  positive  numbers 
are  represented  as  in  the  sign  and  magnitude  format.  A 
negative  number  is  represented  by  complementing  all  of  the 
bits  of  the  positive  number.  In  2's  complement 
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representation,  positive  numbers  are  represented  as  in  the 
sign  and  magnitude  format.  Eut  a  negative  number  is 
represented  by  subtracting  the  magnitude  from  2.0.  The 
choice  among  the  three  formats  depends  on  hardware 
consideration.  The  two's  complement  format  is  widely 
chosen  in  most  digital  systems  because  it  conveniently 
performs  subtraction  using  an  adder.  Another  advantage  of 
using  the  2's  complement  format  is  that  the  correct  total 
sum  will  be  obtained  even  when  partial  sums  overflow  or 
underflow. 

In  fixed-point  arithmetic,  the  result  of  adding  two 
b-bit  numbers  is  still  b  bits.  However,  the  magnitude  of 
the  resulting  sum  can  exceed  unity.  This  phenomenon, 
called  overflow,  is  inherently  related  to  the  limited 
dynamic  range  of  fixed-point  arithmetic.  Scaling  can  be 
performed  to  avoid  undesired  overflow.  The  product  of  two 
b-bit  numbers  results  in  a  2*b-bit  number.  If 
multiplication  is  carried  out  p  times,  the  required 
word-length  for  representing  the  result  is  p-b  bits.  This 
is  clearly  an  unacceptable  situation  for  the  hardware  and 
economy.  To  remedy  this  problem,  truncating  or  rounding 
operations  to  fit  the  results  of  multiplication  into  a 
finite  word-length  register  is  necessary.  The  error  due  to 
truncating  or  rounding  of  p  bits  of  word-length  into  q  bits 
(p  >q)  of  word-length  is  commonly  referred  as  roundoff 

error.  Considerable  attention  has  been  paid  to 
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investigating  the  effect  of  roundoff  error  on  digits] 
filter  implementation  in  the  last  decade  [3-1  to  3-5]. 


Floating-point  arithmetic  is  a  method  for  providing 
automatic  scaling.  An  arbitrary  finite  number  x  can  be 
represented  exactly  using  the  floating  point  representation 

x  =  sign(x)c  •  21  (3-1) 

where  c,  the  mantissa,  is  a  full  precision  binary  number 
such  that  1/2 < c < 1  and  i,  the  exponent ,  is  an  integer. 
The  number  of  bits,  b,  in  a  f 1 owi ng -po i nt  representation 
should  be  divided  into  the  number  of  bits  b^,  for  the 
mantissa  and  the  number  of  bits  b 2  for  the  exponent. 
Although  floating-point  arithmetic  requires  truncating  or 
rounding  operations  in  both  multiplication  end  addition 
[3-6],  it  provides  more  dynamic  range  than  fixed-point 
arithmetic . 

The  comparison  between  fixed-point  and  floating-point 
arithmetic  depends  on  the  input  probability  density 
function,  input  power  spectral  density,  and  filter 
frequency  response  [3-7].  If  the  floating-point  mantissa 
and  fixed-point  word-length  have  the  same  word-length,  then 
floating-point  arithmetic  is  more  advantageous .  Generally, 
when  a  large  dynamic  range  is  required,  floating-point 
arithmetic  generates  less  roundoff  noise  because  it 
provides  automatic  scaling  [3-1].  But  it  should  be  noted 
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here  that  floating-point  arithmetic  is  significantly  more 
complicated  and  expensive  in  hardware  than  fixed-point 
arithmetic.  When  economy  and  speed  are  of  major  concerns, 
fixed-point  arithmetic  is  usually  a  logical  choice. 

The  comparison  between  truncating  and  rounding  depends 
on  whether  fixed-point  or  floating  point  arithmetic  is  used 
and  how  negative  numbers  are  represented.  However, 
experiments  have  shown  that  the  errors  generated  by 
truncation  are  more  severe  than  those  generated  by  rounding 
because  of  a  bias  [3-3].  Truncation  operation  is  not 
commonly  used  in  practical  digital  system. 

The  next  problem  of  concern  is  fraction  or  integer 

multiple  representation  of  numbers.  In  integer  multiple 

-N 

representation,  all  numbers  are  represented  by  2  ,  where  N 

is  an  integer.  Therefore,  the  multiplication  operation 
requires  only  a  shift  operation.  This  shift  operation  will 
increase  computational  speed  and  simplify  the  hardware. 
But  one  can  expect  losses  in  dynamic  range  and  accuracy  in 
arithmetic.  Since  accuracy  is  essential  in  finite 
word-length  arithmetic,  fraction  representation  is  commonly 
chosen  . 

Due  to  all  these  reasons,  attention  will  be  focused  on 
fixed-point  arithmetic  with  the  rounding  operation  and 
fraction  representation.  Another  reason  for  restricting 

our  attention  to  fixed-point  arithmetic  is  that  the 
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overflow  resulting  from  the  limited  dynamic  range  con  be 
avoided  by  proper  scaling  of  the  signal  level. 

3 . 2  Prel iminary  Statement 

Fixed-point  arithmetic  with  finite  word-length  causes 
three  common  error  sources  [3-3]: 

1)  Quantization  of  the  input  signal  into  a  set  of 
discrete  values  causes  inaccuracies. 

2)  Representation  of  the  filter  coefficients  by  a 

finite  word-length  changes  the  filter 

character istics . 

3)  Rounding  or  truncating  of  the  results  of 
arithmetic  operations  within  the  filter  causes 
errors,  known  as  roundoff*  noise  in  the  filter 
output . 

Overflow  can  also  be  a  problem  within  filters.  However, 
the  overflow  problem  can  be  avoided  if  the  signals  are 
properly  scaled.  This  problem  will  be  discussed  later. 

The  first  source  of  error  above,  commonly  referred  to 
as  A/D  noise,  is  inherent  in  any  analog-to-dig ital  (A/D) 


*This  term  is  universally  adopted  whether  rounding  or 
truncation  operation  is  actually  performed. 
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conversion  process,  and  has  been  studied  in  great  depth 
[3-5].  It  is  noted  here  that  the  input  data  array  is 
already  a  quantized  version  in  most  practical  cases.  For 
example,  8-bit  image  data  is  common  in  image  processing. 
Furthermore,  it  shall  be  shown  later  that  the  effect  of 
input  quantization  is  far  less  severe  than  the  effect  of 
roundoff  noise.  Hence,  the  effect  of  A/D  noise  has  been 
excluded  in  this  study. 

The  second  source  of  error  mentioned  above  occurs 
because  the  filter  coefficients,  following  a  design 
procedure,  which  would  normally  use  full  precision,  must  be 
quantized  with  finite  word-length.  This  quantization  of 
the  filter  coefficients  will  alter  the  transfer  function. 
This  error  differs  from  structure  to  structure.  It  is 
advantageous  to  use  a  structure  that  is  insensitive  to 
filter  coefficient  quantization.  In  general,  the  effect  of 
filter  coefficients  in  accuracy  is  most  severe  in  a 
higher-order  filter  when  the  filter  is  realized  in  the 
direct  form  than  when  it  is  realized  in  the  parallel  or 
cascade  form.  As  a  rule,  therefore,  the  parallel  or 
cascade  form  should  be  used  for  higher-order  filters 
whenever  possible  [3-3].  Experimental  results  have  shown 
that  the  amount  of  error  is  not  significant  in  our  case. 
Therefore,  no  particular  emphasis  will  be  made  in  this 
study,  except  in  Chaper  5. 


The  third  source  of  error  mentioned  above  is  of  major 
concern  in  fixed-point  arithmetic,  and  is  the  major  subject 
of  the  next  section.  Roundoff  noise  is  the  most  important 
factor  in  determining  the  complexity  of  hardware  and  speed. 
Large  word-length  will  slow  down  computat ional  speed. 
Furthermore,  the  price  paid  by  increasing  the  word-length 
for  filter  coefficients  is  negligible  compared  to  the  price 
paid  by  increasing  the  word-length  for  reducing  roundoff 
error.  In  addition,  a  limit  cycle  can  occur  in  the 
recursive  realization  of  FIR  filters  [3-9].  However,  a 
limit  cycle  cannot  occur  in  the  nonrecursive  structures. 

3 . 3  Fixed-Point  Arithmetic 

3.3.1  Roundoff  Error 

The  direct  form  of  discrete  convolution  can  be 
characterized  as  a  calculation  of  the  sum  of  products 

N  N 

S  =  ^^a(n)b(n)  =  ^^c(n)  (3-2) 

n=l  n=l 

Let  us  assume  that  a(n)  and  b(n)  are  (b+l)-bit  numbers 
(including  sign  bit)  and  products  are  rounded  to  less  than 
(2b+l)  bits,  but  more  than  (b+1)  bits.  Then,  the  rounded 
products  can  be  written  as 

[c (n) ]  =  c (n)  +  e  (n)  (3-3) 
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The  relation  between  [c(n)]r  and  c(n)  is  shown  in  Fig.  3-1, 
where  [c(n)]  denotes  the  rounded  number  and  e(n) 
represents  the  error  resulting  from  rounding.  In 
fixed-point  arithmetic,  the  error  made  by  rounding  with 
(b^+l)  bits  satisfies  the  inequality 


e  (n)  < 


(3-4  ) 


Thus,  the  resulting  sum  can  be  expressed  as 

N  N 

Sj.  =  [c  (n)  ]  r  =  S  +  ^^e(n)  (3-5) 

n=l  n=l 


Let  us  assume  that  the  resulting  sum  S 
( t>2  +1 )  bits  of  word-length.  Then, 
rounded  to  (fc^+l)  bits  can  be  rewritten 


N 

S2  =  S1  +  v  =  S  +  (n) 

n=l 


where 


_b2  -b2 
2  ^  2 

~2—  *  v^  — 


Therefore,  by  combining  Eqs.  (3-5)  and 


-bf  -b2 

S2'S|  <  N  + 


will  be  stored  into 
the  resulting  sum 
as 


+  v  ( 3~6a) 


(  3  —  6  b ) 


(  3-6 )  ,  we  obtain 


(3-7) 
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The  characteristic  of  the  roundoff  noise  at  the  output 
depends  on  the  location  where  rounding  is  performed.  There 
are  two  possible  locations  for  rounding.  First,  if  all 
mul t ipl ica t ions  are  performed  with  full  precision,  rounding 
is  performed  only  after  summation.  Then,  from  Eq.  (3-5), 
e(n)=  0  for  n  =  1,2, ...N  so  that 

2~b2 

|S2-S  |  <  —  (3-8) 

multiplication  products  are  rounded  for  storage 
addition,  v  =  0  and 

2~bl 

|S2-s|<£_n  (3-9) 

Unfortunately,  all  of  the  bounds  derived  are  for  worst 
cases,  and  thus,  are  of  little  practical  usefulness.  In 
the  following  discussion,  we  will  derive  more  useful 
bounds . 

A  less  conservative  estimate  of  the  noise  caused  by 
rounding  can  be  obtained  by  a  statistical  approach  (3-3]. 
It  should  be  noted  here  that  a  precise  analysis  of  roundoff 
noise  is  generally  complicated,  and  not  reauired  in 
practical  appl icat ions .  The  purpose  of  error  analysis  is 
to  determine  word-length  within  a  filter  to  satisfy  some 
specification  with  reasonable  tolerance.  Furthermore,  a 


If  all 
before 
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final  decision  concerning  word-]ength  is  insensitive  to 
inaccuracies  in  the  error  analysis.  Thus,  an  analysis 
correct  to  within  30  %  to  40  %  is  often  acceptable  (3-6). 

The  statistical  approach  considers  the  errors 
introduced  by  rounding  to  be  small  random  quantities.  This 


viewpoint 

simpl  if  ies 

the 

analysis  and 

enables 

useful 

theoretical 

resul ts 

to 

be  derived . 

Many 

computer 

simulat ions 

results 

have 

verifed  the 

val id  ity 

of  the 

statistical 

appr oach 

[3-3 

to  3-5]  .  It 

has  been 

claimed 

that  the  statistical  approach  tends  to  be  more  accurate 
when  the  number  of  quantization  levels  is  not  too  small 
[3-3]  . 

Three  common  assumptions  are  made  concerning  the 
effect  of  rounding  [3-3].  They  are: 

1)  The  error  sequence  e(n)  is  a  white-noise  sequence. 

2)  The  probability  distribution  of  the  error  sequence 
e(n)  is  uniform  over  the  range  of  quantization 
intervals . 

3)  The  error  sequence  e(n)  is  uncorrelated  with  the 
input  and  itself. 

The  uncorrelatedness  assumption  is  particularly 

attractive  because  the  total  error  due  to  rounding  is  the 

sum  of  each  rounding  error.  There  are  some  controversies 
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OVer  '>»  validity  „f  th 

« -  tBm  is  c  n;;;;;ptio"s  -  - 

r au  :  -» - 

CSSeS'  the  roundoff  noise  ,  ^  lnval‘<>-  in  such 

iS  ~  ^  uncorrelated  WiVT^'  ^ 

;;;;;tions  —  *  -  —  for  most  efi;;put-  -  — 

signals  of  z  filters  uit-h 

reasonable  ampiitude  '"Put 

uncor relatedness  is  not  »««t.  If 

"’°re  and  the  '  "  ^  be 

‘“'“CUl«  !"P-  signal  or  class  f^  ^  *"’”*”*  «  *h, 

Clasa  ot  input  signal s  f3_,, 

8“-  on  the  discussion  above  p1q  , 

“  r'Pl3Ced  by  “  «««..  roundoff  nof  °P*r«'<>» 

“e  aSSU",s!d  ^at  all  multipi  ication  ^  this  m°d*1  ■ 

o«ctly,  snd  round.ng  Js  per  Pr°dU“s  ate  represented 

I-.,  at  the  fiifet  ^  are 

source  ls  present  Jn  ^  f  .  ,  ^  ’  Then  ""'i  one  noise 

°“tput-  '  and  “  soper imposes  on  the 

There  3re  on  ?  i 

6  20  3xl  SGK  filter*  ■ 

stage ,  n  ln  each  svd  o 

C  columns  and  W  expansion 

V  r  O  WS  ,  r  _  i. 

Let  “S  define  a 


two-dimensio„ai  3x3  fllt,r  t  ^  Ut  ™  *«„.  . 

col“»  or  ro»  direction.  '  depe"d1"^  on  the 

SVD  expansion.  Thus,  Ut  ^  3  «»  J-th 
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t .  .  ( l ,m) 

3  > 1 


0' 

0 

0 

01 

r. 

0 


for  column  convolution 


(3-10) 


for  row  convolution 


for  i=l , 2 , . . . , 2Q .  Figure  3-3  shows  the  roundoff  noise 
model  for  the  SVD/SGK  convolution  filter.  The  mean  and 
variance  of  the  error  sequence  e(n)  can  be  shown  to  be 


m  =  0 
e 


(3-11) 


>-2b 

TT~ 


We  assume  that  the  rounding  is  performed  with  (b+1)  bits 
word-length.  in  this  model,  a  given  error  sequence  e(n)  is 
filtered  by  succeeding  sections,  so  that  the  output  noise 
variance  will  depend  upon  the  ordering  of  the  second  order 
SGK  filters. 


Let  us  define  g.  .(£,ra)  to  be  the  impulse  response  from 

l,i 


the  noise  source  e^(n)  to  the  output.  Thus, 


gj,i(^,m)=tj  ,i+l<?',rn)®®tj  ,i+2(<i'm)99'  '  j  ,2Q(!l'ro)  (3“12) 


The  mean  and  variance  of  the  roundoff  noise  are  then  given 
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by 


m 


e 


i 


0 


zv 


a  . ( £ ,m) 
3  /  ! 


m 


(3-1 ?) 


and  the  total  noise  variance  is  the  sum  of  each  noise 
variance  of  the  3x1  SGK  filter.  Therefore, 


2Q 

i=l  £  m 

If  an  impulse  response  H  is  approximated  by  K  singular 
values,  then  the  total  output  noise  variance  due  to 
rounding  is 


(3-15) 


If  the  two-dimensional  impulse  response  g.  .  (£,m)  consists 

Jf  1 

of  N  ^  SGK  filters  for  the  columns  and  (2C-i-Nj_)  SGK  filters 

for  the  rows,  then  g.  .(£,m)  can  be  rewritten  as 

3  > 1 


g  . (£,m) 
3  '  A 


gi  t(£)gi  i{m) 
J  t  *  J  /  X 


(3-16) 
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In  vector  form,  Eq.  (3-16)  is  equivalent  to 

c  ,  r  .  T 
g,  i  =  i(g-i  J 

~J  /  X  J  9  J-  J  J1 


(3-17) 


where  gc  .and  .  are  one-dimensional  impulse  responses 
obtained  by  convolving  SGK  filters  for  the  columns,  and 
( 2Q-i-N, )  SGK  filters  for  the  rows,  respectively.  If  g 
consists  of  only  SGK  filters  for  the  columns  or  the  rows , 
then  gc  and  gc  should  be 

j/i  j/i 


(3-1 8a) 


(£4  i)T  =  to  1  0] 

“J  >  x 


(3-1 8b) 


Note  that 


Am  l  m 


(3-19) 


Substituting  Eq.  (3-19)  into  Eq.  (3-20),  we  obtain 


K  2  Q 


’total"  TT  E{  £  [E  lgj,iU),2£,gj,i(m)  |2]|  (3‘ 

j=l  i=l  i  m 


64 


Equation  (3-20)  is  a  theoretical  formula  predicting 
roundoff  noise  with  (fc+1)  bits  word-length.  Tts  validity 
will  be  demonstrated  in  Chapter  5. 

3.3.2  A/D  Error 

Next  an  attempt  has  been  made  to  show  that  the  input 
A/D  noise  is  negligible  compared  to  the  roundoff  noise. 
Again,  the  statistical  model  is  chosen,  and  the  input 
quantization  is  considered  as  an  injection  of  additive 
noise  to  the  input.  The  noise  cuantitites  are  uniformly 
distributed  over  one  quantization  interval  and 
statistically  independent.  Since  the  first  place  where 
quantization  of  the  input  signal  may  take  place  is  at  the 
A/D  converter,  the  A/D  noise  effect  is  independent  of  the 
structure  we  used  to  realize  the  filter.  Figure  3-4 
describes  the  statistical  model  for  A/D  noise. 

If  the  quantizer  has  a  word-length  of  (t  +  1)  bits,  then 
the  input  to  the  actual  filter  is  x  (  £,m)  £,m)  ,  where 

e  (  9  ,m)  is  the  quantization  error,  bounded  by 

n/  U 

2  —  t  o  —  t 

- 0—  <  e  (  £,m)  <  t _ .  Let  us  define  the  output  error 

1  A/D  2 

array,  E ( 9  ,m)  ,  as 

E(2,m)  =  H ( i ,m)®@eA^D ( l ,m)  (3-21) 

Since  the  filter  is  linear,  it  can  be  shown  that  Fltfm)  has 
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zero  mean  ana  variance  given  by 


2 

°A/D 


'l  EE 


H  (  '«  ,  m ) 


m 


where 


a 


2 

t 


( 3-2  2a ' 


( 3-2  2b) 


It  is  noted  that  the  filter  has  been  normalized,  so  that 


H ( ? ,m) 


1 


(3-23) 


Such  a  normalized  filter  will  not  change  image  contrast 
between  input  and  output.  Therefore,  it  is  obvious  that 


EE  |HU,m)  |2<  EE  HU,m)  =  1  (3-24) 

l  m  i  m 

Using  Eqs .  (3-17)  and  (3-19),  and  assuming  that  the 
quantizer  has  the  same  word-length  as  a  multiplier,  it  can 
be  shown  that 

°A/D  S  °e  «3-25' 

It  is  shown  that  the  A/D  noise  is  smaller  than  or  equal  to 
that  of  roundoff  noise.  in  general,  A/D  noise  is 
negligible  compered  to  roundoff  noise. 


It  is  necessary  to  remark  on  the  effect  of  filter 
coefficient  quantization.  Although  zero  location  and 
frequency  response  sensitivities  to  coefficient  changes  can 
readily  be  obtained,  no  general  statistical  analysis  of  the 
type  given  in  Section  3-3  has  been  obtained  to  describe  the 
cascade  form  of  FIR  filters.  Herrmann  and  Schuessler  [3-91 
worked  on  this  problem  only  experimentally,  not 
theoretical  1 y . 

3 . 4  Conclusion 

The  accuracy  of  a  digital  filter  is  limited  by  the 
finite  word-length  used  in  its  implementation.  When  a 
digital  filter  is  implemented  with  special-purpose 
hardware,  one  is  usually  interested  in  determining  the 
minimum  word-length  needed  for  a  specified  performance 
accuracy.  Also,  word-length  is  an  important  factor  in 
determining  the  complexity  of  hardware  and  speed.  Thus,  it 
is  very  important  to  understand  the  effect  of  quantization. 

In  this  chapter,  attempts  have  been  made  to  analyze 

relevant  effects  of  using  fixed-point  arithmetic  for 

SVD/SGK  convolution  from  a  statistical  viewpoint.  Cur 

consideration  of  finite  word-length  effect  began  with  a 

discussion  of  the  various  methods  of  number  representation 

that  are  commonly  used  in  digital  system.  The  following 

discussion  focused  on  three  common  source  of  errors  caused 

by  implementation  with  finite  word-length.  Then,  we  showed 
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how  a  statistical  analysis  can  estimate  the  effects  of 
quantization  in  SVD/SCK  convolution.  Statistical  methods 
were  shown  to  be  very  efficient.  in  systems  with 
non-determ inist ic  signals.  It  was  also  shown  that  roundoff 
noise  is  of  major  concern  in  digital  implementation,  and  a 
theoretical  formula  to  predict  total  roundoff  noise 
variance  of  SVD/SGK  convolution  was  derived.  The  A/D  noise 
was  shown  to  be  negligible  in  our  case.  Finally,  the 
dependence  of  the  roundoff  noise  on  section  ordering  was 
demonstrated.  The  discussion  of  section  ordering  and  a 
dynamic  range  consideration  in  the  fixed-point  arithmetic 
is  the  subject  of  the  next  chapter. 
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CHAPTER  4 


SCALING  AND  SECTICN  ORDERING 


4.1  Introduction 


In  the 

pr ev ious 

chapter  ,  a 

theoretical 

formula 

for 

predicting 

roundoff 

noise  variance  was  derived. 

Cne 

important 

constraint  should 

be  imposed 

on 

the 

implementat 

ion  with 

fixed-point 

arithmetic . 

There 

is  a 

finite  dynamic  range 

of  fixed-point  arithmetic. 

To  ensure 

that  the  final  output  be  correct,  overflow  at  the  output  of 
any  second-order  SGK  filter  must  be  avoided.  if  the  output 
of  each  section  (SGK  filter)  exceeds  the  finite  dynamic 
range  of  the  filter,  undesired  signal  distortion  will  be 
introduced  to  the  output.  For  example,  given  the  dynamic 
range  of  (-1.0,  1.0),  adding  two  numbers  may  result  in  a 
number  that  is  not  within  the  given  range.  Truncating  or 
rounding  operations  that  assign  the  limit  value  to  the 
result  (say  -1.0,  or  1.0)  introduce  an  error.  This  problem 
directs  attention  to  the  need  for  a  scaling  procedure  for 
the  filter  parameters  of  each  SVD/SGK  section  in  order  to 
prevent  overflow. 

Another  issue,  section  ordering,  is  also  important  to 
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minimize  roundoff  noise.  As  seen  in  Fq.  (3-15),  the  total 

roundoff  noise  has  a  strong  dependence  on  section  ordering, 

i.e.,  g  (>°  ,m)  will  be  different  if  the  ordering  is 
j/i 

different.  For  example,  Schussler  [4-]]  has  demonstrated 

that  a  FIR  filter  with  a  length  of  33,  ordered  one  way, 
2  2 

produces  a  =  2.40  ,  where  Q  is  the  quantization  step, 

2  g  7 

while  ordered  another  way  yields  o  =  1.5x10  Q' .  Although 
this  experiment  demonstrated  two  extreme  cases,  it  clearly 
shows  the  importance  of  section  ordering  in  cascade  form 
FIR  filters.  Since  the  respective  difference  is  so  large, 
determining  the  minimum  roundoff  noise  ordering  is 
essential  . 

Unfortunately,  attempts  to  find  optimal  ordering 
become  impractical  since,  given  M  sections,  there  are  M! 
possible  orderings.  Even  for  s  moderate  value  of  M,  say 
M  =  7,  searching  5040  possible  orderings  is  very 
time-consuming.  Furthermore,  due  to  the  analytical 
complexity  of  Fq.  (3-15),  no  analytical  approach  to  finding 
an  optimal  ordering  seems  possible. 

Chan  and  Rabiner  [4-2]  investigated  the  section 
ordering  problem  of  one-dimensional,  cascade  form  FIF 
filters  quite  intensively  and  reported  their  results,  based 
on  the  experiment,  as  follows: 

i)  Most  orderings  have  very  low  noise  compared,  to  the 
maximum  possible  value.  More  specific,  for  a  FTP 
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lowpass  filter  with  length  31,  they  showed  that 
approximately  two-thirds  of  the  orderings  have 
noise  variance  less  than  4%  of  the  maximum,  of 
which  nine-tenth  have  noise  variance  less  than  14% 
of  the  maximum. 

2)  There  is  a  large  gap  between  small  and  large  noise 
variance  distribution,  and  the  noise  values  within 
the  gap  are  produced  by  very  few  orderings. 


Their  conclusions  are  encouraging.  Since  the  larqe 
majority  of  possible  orderings  are  very  close  to  the 
minimum  noise  variance  ordering,  finding  a  suboptimal 
ordering  is  possible  with  reasonable  computations.  Instead 
of  finding  a  time-consuming  optimal  ordering,  it  may  be  far 
more  practical  to  use  a  suboptimal  ordering  method  that  can 
rapidly  determine  an  ordering  that  is  close  to  the  optimal. 
Furthermore,  the  reduction  in  roundoff  noise  gained  by 
finding  the  optimum  ordering  is  very  small,  compared  to  a 
good  suboptimal  ordering. 

Based  upon  their  experiments,  Chan  and  Fabiner 
proposed  a  simple  one-dimensional  ordering  algorithm  [4-21, 
which  has  proven  to  be  very  efficient  in  minimizing 
roundoff  noise  variance. 

The  purpose  of  this  chapter  is  to  discuss  seeling 
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procedures  and  ordering  algorithms  for  SVD/SGK  convolution. 

Two  existing  scaling  methods,  sum  and  L  -norm  scaling,  are 

P 

discussed,  and  applications  to  SVD/SGK  convolution  are 
given  in  Section  4-2.  A  brief  review  of  the  Chan  and 
Rabiner  ordering  algorithm  and  its  generalization  to 
SVD/SGK  convolution  are  described  in  Section  4-3. 

4 . 2  Seal ing  Procedure 

Scaling  is  important  because  the  computational  dynamic 

range  sets  a  practical  limit  to  the  maximum  value  of  signal 

levels  representable  in  the  filter.  The  theoretical  basis 

for  the  scaling  procedure  chosen  here  is  Jackson's  work 

[4-3],  commonly  referred  to  as  sum  scaling.  To  formulate 

the  required  overflow  constraints,  let  us  assume  that  an 

input  signal  x(n,m)  is  bounded  in  magnitude  by  1.0. 

•  '  •  *  •*  . . .  .  . 

We  shall  consider  a  scaling  procedure  in  which  a 

( 20+1 ) x ( 2Q+1 )  FIR  filter  is  implemented  by  SVD/SGK 

convolution.  There  are  2Q  3x1  SGK  filters  for  the  columns 

and  rows  in  each  SVD  expansion  stage.  To  simplify 

notation,  only  one  SVD  expansion  stage  will  be  considered. 

Therefore,  the  subscript  j  will  be  dropped. 

We  will  define  f^(2  ,m)  to  be  the  impulse  response  from 
the  input  to  the  i-th  section.  The  z-transform  of  f^(2  ,m) 


can  be  written  as 


(4-1  ) 


Fi(zi'z2}  =£i;fiu'm)ziS"m 

i  m 


n  Tp(zi'z2^ 

P=1  ^ 


where  Tp(z-L,z2)  is  the  z-transform  of  tp  (  5.,m)  ,  as  defined 
in  Eq.  (3-10).  Let  be  the  scaling  factor  for  the  i-th 

I 

section  and.  ( z  ,z^ )  be  a  scaled  z-transform  of  T\(z^,z2). 

Then 

t!_(z1,z2)  =  SiTi(z1,z2)  (4-2) 

and  the  scaled  transfer  function  from  the  input  to  the  i-th 
section  is 


Fi<zi'*2)  fi  U,m)  z^'z 


-£_-m 
2 


m 


( 4  —  3  e ) 


or 


Fi(zl'z2) 


i  ,  li 

n  n  s .  n  r^r-, ( z i • z ? ) 

p=i  p  p=i  1  p=i  F 


( 4  —  3  b ) 


Letting  v^(?.,m)  be  the  output  at  the  i-th  section,  v^(jj,,m) 
is  obtained  by  convolving  the  input  array  x(2,m)  with  the 

I 

impulse  response  f^(jl,m).  Thus 

SL  m 

vi(!l,m)  =  x(p,q)  f  .  U-p+l,m-q+l) 

p  q 


(4-4) 


where  g'(u,v)  and  i,(u,v)  are  the  Fourier  transform  of 

I 

and  x(£,m)  respectively.  Fere  we  assume  that  the 
input  x(2,m)  is  a  deterministic  signal. 


for 


1 

P 


+ 


1 


and  p,q>l.  Therefore,  with  L^-norm  scaling,  the  reauired 
condition  and  preventing  overflows  is  satisfied  by 


1  v  i  (  £  ,  m )  |  <  ||f’H  p 


(4-1 ?) 


Based  on  Eq.  (4-13),  a  sufficient  condition  for  scaling  can 
be  given  if  one  has  knowledge  of  the  Lp-norm  of  the  input 
signal.  One  particularly  interesting  case  is  p  =  and 
q  =  1.  In  this  case,  the  input  signal  should  be  bounded  by 


;k/7 


|  X  (u , v) | dudv  <  1.0 


-IT  -IT 


(4-14) 


Then,  the  necessary  and  sufficient  condition  on  the  scale 
factor  to  guarantee  that  the  output  of  each  section  is* 
bounded  in  magnitude  by  1.0  is  that 


II 


t 
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-7T<U<TT 
—  TT<  VfCTT 


(u,v) | 


< 


1.0 


(4-15) 


which  is  equivalent  to 


n  s .  < 

j=i  : 


MAX  , 

-TT<U<TT  |  '3  .  (u  ,  v) 
-TT<V<TT  1 


(4-16) 
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The  scaling  procedure  of  Eg.  (4-16),  satisfied  with 
equality,  is  known  as  peak  scaling  [4-3]. 

The  two  scaling  procedures  discussed  above  are 
summarized  as  follows.  If  the  input  signal  is  bounded  by 
|  x(£,m)|  <  1.0,  then  each  scale  factor  can  be  compi  ed 
according  to  the  following  procedure: 


1.  Compute 


a?  =^^|fi(£,m)|  for  i=l,2,...,2Q 
a  m 


(4-1 7  a) 


2.  Then 


S.  = 

l 


Vi 


for  i=l 


for  i=2 , 3 , . . . , 2Q 


If  input  signal  is  bounded  such  that 


TT  TT 

I1*'"''' 


)  |  dudv  <  1.0 


-IT  -TT 


(4-1  7b) 


(4-1 8a) 


then  each  scale  factor  is  computed  as  follows: 
1.  Compute 
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for  i=l , 2 , . . .  ,  2Q 


{ 4  —  1 8  b ) 


y 


MAX 

|j.  (u,v) 

-  7T<V  <TT  1 


2.  Then 


for  i=l 


for  i  =  2 , 3  ,  .  .  .  ,  2Q 


(4-1 8c) 


But,  it  should  be  noted  here  that  Eq.  (4-16)  cannot  be  used 
in  the  case  of  a  random  input  signal,  because  if  the  signal 
is  random,  its  Fourier  transform  does  not  exist..  Instead 
of  using  the  Fourier  transform,  the  equivalent  condition 
can  be  obtained  with  the  appropriate  power  spectral  density 
and  autocorrelation  function  [4-3). 


Experimental  results  indicate  that  the  two  scaling 

•  •  •  • 

procedures  yield  noise  variances  that  are  similar  (4-21. 

In  general,  sum  scaling  is  much  simpler  to  perform  than 

peak  scaling  in  FIR  filter  cases.  With  peak  scaling,  one 

must  find  the  maxima  of  the  |3j_(u,v)  |  for  all  i.  Even 

using  the  FFT  algorithm  will  require  more  computations  than 

finain,  EE  |  f  i  (  2,  ,m)  I  for  all  i.  It  has  been  claimed 

l  m 

that  sum  scaling  is  too  conservative  to  be  used  in  TIP 
filter  cases  [4-4].  But  this  is  not  of  major  concern  in 
the  case  of  a  FIR  filter.  Therefore,  in  order  to  save 
computation  time,  we  shall  focus  on  sum  scaling  throughout 
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this  study. 


The  sum  scaling  of  Fa.  (4-7)  reauires  computation  of 
the  two-dimensional  impulse  response  f^(i,m)  for  all  i. 
Since  each  SVD-expanded  matrix  H  of  Fa.  (2-]  7)  is 
separable,  further  simplification  is  possible  for  FVD/SCK 
convolution.  Note  that  each  separable  matrix  F  is  an 
outer  product  of  one-dimensi.onal  column  and  row  convolution 


operators  c  .  and  r  .  . 

~3 


Instead  of  applying  the  sum  seal  ing 


by  computing  f^(H,m),  the  same  result  will  be  obtained  by 


applying  the  sum  scaling  to  c.  and  r.  independently.  The 

"I  ~J 

following  Lemma  will  generalize  the  above  argument. 


Lemma ;  If  a  two-dimensional  separable  impulse  response 
matrix  H  is  realized  in  the  one-dimensional  cascade 
form,  sum  scaling  can  be  applied  to  the 
one-dimensional  convolution  operators  c  and  r 
independent! y . 


Proof :  In  the  SVD/SGK  convolution  system,  there  are  0 

3x1  SGK  filters  for  the  columns  and  rows  of  the 
input  image.  Given  a  certain  ordering,  there  are  N 
3x1  SGK  filters  for  the  columns  and  (i-N^)  3  x  1  SGK 
filt  ers  for  the  rows  from  input  to  the  i-th  section. 
Let  us  assume  the  i-th  section  is  a  filter  for  the 
column.  From  the  sum  scaling  of  Fq.  (4-17),  we  have 
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(4-19) 


which  is  equivalent  to  one-dimensional  sum  scaling. 

By  the  Lemma,  two-dimensional  sum  scaling  is  shown  to 
be  equivalent  to  two  one-dimensional  sum  scaling 
operations.  The  same  Lemma  can  be  applied  to  peak  scaling. 
Since  fi(2,m)  is  separable,  then 

3i(u,v)  =  ^  (u)  J  (v)  (4-22) 
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Therefore 


"  MAX 

-TT<U<TT  |  3K  (u,  V)  | 
-TI<V<7T 


MAX 

j-Tt<U<TT 


MAX 

-TT<V<7T 


(4-23) 


4 . 3  Section  Ordering 


The  next  step,  given  the  scaling  procedure  chosen,  is 
to  choose  an  ordering  for  the  sections  to  minimize  the 
total  roundoff  noise.  As  an  approach  to  determine  an 
ordering  algorithm  for  SVD/SGK  convolution,  a 
one-dimensional  ordering  algorithm  for  the  cascade  form  FIR 
filter  will  be  introduced.  If  a  FIR  filter  of  size  (20+1) 
is  realized  in  cascade  form,  there  are  0  sections  of 
second-order  filters  and  C!  possible  orderings.  If  we 
define  b^(k)  for  i=l,2,...,C  to  be  the  impulse  response 
from  the  (i+l)-st  section  to  the  output,  the  total  roundoff 
noise  variance  can  be  shown  to  be  [4-2] 


•2b 


a  = 


12 


Q 

ECE 


|bi(k) 


i=l 


(4-24) 


Here  we  assume  that  the  rounding  is 
the  products  are  represented  in 
ordering  will  minimize  the  total 
Based  on  the  Chan  and  Fabiner 
algorithm  is  summarized  as  follows 


performed  only  after 
full  accuracy.  The  best 
output  noise  variance, 
experiment,  the  proposed 
[4-2]  : 
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i-th 


section 


Beginning  with  i  =  Q,  assign  the 

together  with  the  section  already  assigned,  that 

yields  the  smallest  possible  value  for  b^_1  (k)|^. 

k 

This  algorithm  is  suboptimal  since  the  algorithm 
minimizes  the  output  noise  variance  from  individual 
sections  instead  of  minimizing  the  sum  of  the  output  noise 
variance.  However,  'in  all  cases  tested,  the  algorithm  has 
proved  to  yield  section  ordering  very  close  to  the  optimum 
ordering  because  a  large  majority  of  possible  orderings 
yield  small  output  noise  variance,  as  discussed  before. 

In  SVD/SGK  convolution,  there  are  a  total  of  20  3x1 

SGK  filters.  Searching  (2Q) !  possible  ordering  is  an 
enormous  task.  But  we  shall  see,  based  on  the  existing 
theory,  the  generalization  of  a  one-dimensional  ordering 
algorithm  for  SVD/SGK  convolution  is  possible,  and  the 
proposed  algorithm  will  prove  to  be  efficient  and  simple. 


Let  us  rewrite  the  output  noise  variance  formula  for 
SVD/SGK  convolution,  as  derived  in  Eqs.  (3-15)  and  (3-20), 
as 


where  g  (  £,m)  ,  g.c.(£),  and  g.r.(m)  are  already  defined  in 
Ji  1  J  /  1  3  /  1 

Section  3-3.  Again,  only  one  SVD  expansion  term  will  be 
considered;  therefore,  the  subscript  j  will  be  dropped. 


Using  Eg.  (4-25),  it  is  quite  simple  to  extend  the 
Chan  and  Rabiner  ordering  algorithm  to  SVD/FGK  convolution. 
Eut  the  significance  of  using  Eq.  (4-26)  to  search  for  an 
ordering  algorithm  for  SVD/SGK  convolution  is  that  the 
ordering  problem  can  be  treated  as  solving  two 
one-dimensional  ordering  problems.  Since  ^  |g^(£)|2  and 


£  K 

m  x 


(m) 


are  positive  numbers,  the  following  is 


satisfied: 

min 
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Rather  than  minimizing  7.7.  |g  (£,m) 

Tm  i 


(4-27) 


an  equivalent 


?!•? 


(£>| 


and 


condition  can  be  obtained  by  minimizing 

1^  separately.  Thus,  minimizng  |9^(£)|”  and 

r  2  £^ 

2.|9i(m)|  is  equivalent  to  two  one-dimensional  ordering 

nr 

problems.  After  ordering  column  and  row  operators 
independently,  the  remaining  step  is  to  decide  whether  the 
SGK  filter  on  the  column  or  on  the  row  should  be  assigned 
at  the  i-th  section. 


To  show  the  rationale  for  the  algorithm 
mathematically,  assume  that  { 2Q— i )  SGK  column  operators  and 
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col umn 


row  operators  have  been  already  assigned  with  F 
operators  and  (20-N-^-i)  row  operators.  Now,  we  want  to 
select  the. i-th  section  of  the  5YD/SGK  convolution.  To 
simplify  the  discussion,  let  us  define 

ai  =  H  lgi+l  a)  I2  (4-28) 

£ 

Bi  lgi+l(m)  I2  (4-29) 

m 

If  we  had  assigned  the  next  filter  on  the  columns  to  the 
i-th  section,  then  the  output  noise  variance  would  be 
proportional  to 

Ei  *  Bi  E|giU)|2  (4-30) 

£ 

If  we  had  assigned  the  next  filter  on  the  rows  to  the  i-th 
section,  the  resulting  output  noise  variance  would  be 
proportional  to 

Ei  =  ou  £|g*(m)  l2  (4-31) 

m 

By  comparing  E ^  and  E  ^  of  Eqs.  (4-30)  and  (4-31),  one  can 
easily  decide  whether  the  filter  on  the  columns  or  the 
filter  on  the  rows  should  be  assigned  to  the  i-th  section. 

Since  0^,8^  for  i=l,2,...C  can  be  obtained  as  a 
result  of  one-dimensional  orderings  for  the  column  and  row 
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operators,  this  procedure  is  far  simpler  than  using 
Eq.  (4-25). 

In  brief,  the  proposed  ordering  algorithm  is 
summarized  as  follows: 

1.  Find  a  one-dimensional  ordering  to  the  column  and 
the  row  operators  by  using  the  Chan  and  Babiner 
algorithm  and  store  c^, 0^  for  i=l,2,...0. 

.  .  .  c  r 

2.  Beginning  with  i=2Q,  compare  and  E.  given  by 

Eqs.  (4-30)  and  (4-31),  respectively,  to  decide 
whether  the  filter  on  the  rows  or  the  filter  on 
the  columns  should  be  assigned  to  the  i-th 

section . 

This  proposed  algorithm  is  also  suboptimal  in  minimizing 

7;  £|g.U,m)|2  rather  than  minimizing  X/  a  •  >  where 
l  m  1  i  1 

°i  *  ££  |9,  <d,m)  j2. 
dm 

4.4  Conclusion 

In  addition  to  the  effect  of  finite  word-length 
discussed  in  Chapter  3,  the  problems  of  overflow  and 
section  ordering  to  minimize  the  total  roundoff  noise  are 
of  great  importance  when  a  digital  filter  is  realized  in 
cascade  form.  To  prevent  overflow,  the  filter  parameters 
and  input  signals  must  be  scaled  so  that  no  overflow  occurs 

following  addition.  Proper  ordering  must  also  be  found  for 
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a  filter  in  cascade  form  because  the  output  roundoff  noise 
has  a  strong  dependence  on  the  way  it  is  ordered. 

Following  the  discussion  of  two  different  scaling 
methods,  sum  and  Lp-norm  scalings,  sum  scaling  was  chosen 
because  sum  scaling  is  simple  and  easy  to  employ.  A 
detailed  sum  scaling  procedure  for  SVD/SGK  convolution  was 
presented.  Because  separable  matrices  result  from  the  SVD 
expansion  of  a  nonseparable  impulse  response  matrix,  the 
two-dimensional  scaling  problem  turned  out  to  be  two 
one-dimensional  scaling  problems.  The  proof  was  given  in  a 
Lemma . 

Next,  the  section  ordering  problem  was  considered. 
Extending  the  existing  one-dimensional  suboptimal  ordering 
algorithm  proposed  by  Chan  and  Rabiner  [4-4],  a  generalized 
two-dimensional  suboptimal  ordering  algorithm  for  SVD/FGK 
convolution  was  proposed.  Because  it  is  actually 
equivalent  to  two  one-dimensional  ordering  problems,  the 
proposed  ordering  algorithm  is  very  simple  and  fast  to 
compute.  The  experimental  results  based  on  the  proposed 
ordering  algorithm,  which  is  shown  to  be  very  efficient, 
will  be  discussed  in  Chapter  5. 


CHAPTER  5 


EXPERIMENTAL  RESULTS  OF  SVD/SGK  CONVOLUTION 
USING  FIXED-POINT  ARITHMETIC 

5. 1  Introduction 

In  this  chapter,  computer  simulation  experimental 
results  for  SVD/SGK  convolution  are  presented.  Two 
prototype  filters  with  linear  phase  have  been  chosen  for 
the  experiments.  One  is  a  lowpass  filter,  the  other,  a 
bandpass  filter.  The  sizes  of  the  filters  are  15x15  and 
11x11,  respectively.  Perspective  views  of  the  frequency 
response  of  the  prototype  filters  are  given  in  figures  5-1 
and  5-2,  respectively.  Figure  5-3  shows  the  SVD 
approximation  errors  for  the  expansion  of  the  prototype 
filters.  It  is  observed  that  the  NMSE  decreases  very 
rapidly  in  both  cases.  In  the  case  of  the  lowpass  filter, 
the  SVD  approximation  error  with  3-stage  expansion  is 
0.5336  %.  In  the  case  of  the  bandpass  filter,  the  SVD 
approximation  error  for  a  4-stage  expansion  is  0.7825  %. 
Numerical  and  photographical  results  related  to  the  outputs 
of  this  SVD/SGK  convolution  when  the  inputs  are  random 
numbers  and  real  image  are  presented  in  this  chapter. 
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Figure  5-1. 


Perspective  view  of  the  frequency 
response  of  the  prototype  bandpas 
filter 


NO.  OF  SVD  TERMS 


Figure  5-3.  NMSE  versus  number  of  singular  value 


5 . 2  Fixed-Point  Ar i thmet ic  Exper iment al  Results 


Our  experiments  were  made  in  the  following  framework. 
We  shall  use  M  to  denote  the  word-length  for  representing 
filter  coefficients  and  N  to  denote  the  word-length  for 
storing  intermediate  results.  Furthermore,  we  shall  adopt 
the  policy  that  all  signal  levels  that  are  representable  by 
given  word-length  be  constrained  within  the  range  of 
(  -1,1).  A  multiplier  with  input  signal  level  greater  than 
unity  may  need  to  be  followed  by  extra  accumulators  and 
extra  wide  adders.  Hence,  more  hardware  is  reauired.  The 
number  of  rounding  operations  within  3x1  SGK  filters  is 
assumed  to  be  one.  In  other  words,  since  the  typical 
operation  performed  in  convolution  is  a  sum  of  products,  we 
assume  here  that  the  rounding  operation  is  performed  only 
after  the  products  have  been  summed  with  full  precision. 
In  addition  to  this,  the  cascade  form  of  the  SVD/SGK 
convolution  requires  a  proper  section  ordering.  The 
suboptimal  ordering  algorithm  discussed  in  Chapter  4  was 
adopted  to  minimize  roundoff  noise.  Because  of  the 
quadrilateral  symmetry  of  the  prototype  filters  used,  the 
one-dimensional  column  and  row  convolution  operators 
obtained  from  the  SVD  expansion  of  H  were  identical.  Thus, 
their  cascade  forms  were  identical.  The  ordering  algorithm 
ended  with  a  perfect  interlace  scheme;  each  filter  for  the 
rows  convolution  was  followed  by  a  filter  for  the  columns 
convolution  and  vice  versa.  But  it  can  be  proved  that  this 
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result,  although  often  true,  cannot  be  generalized  to  all 
esses  of  quadrilateral  symmetrical  filters. 

Then,  after  the  ordering  procedure,  sum  scaling  was 
applied  to  the  filter  coefficients  so  that  overflow  will 
not  occur  within  filters.  Unfortunately,  large  differences 
in  magnitude  among  the  coefficients  causes  the  scaled 
filter  coefficients  to  exceed  the  given  dynamic  range  of 
word-length.  In  this  case,  the  filter  coefficients  were 
further  divided  by  their  maximum  coefficient  to  insure  that 
the  scaled  filter  coefficients  lie  within  the  given  dynamic 
range  of  word-length  for  the  filter  coefficents. 

5.2.1  Roundoff  Error 

To  confirm  the  validity  of  the  noise  formula  of 
Eq.  (3-15),  derived  in  Chapter  3,  a  uniform  density  random 
number  array  of  46x46  pixels  has  been  generated  as  an 
input.  The  statistical  approach  used  to  analyze  roundoff 
noise  in  Chapter  3  is  not  practical  if  the  input  is 
deterministic.  For  this  analysis,  an  image  array  has  been 
modeled  as  a  Markov  process  with  an  adjacent  pixel 
correlation  coefficient  along  lines  of  0.95  [5-1]. 
Furthermore,  it  has  been  assumed  that  the  maximum  signal 
magnitude  of  the  input  array  is  unity,  so  that  all  signals 
are  represented  by  given  dynamic  ranges  of  word-length. 

If  the  filter  size  is  15x15,  then  the  output  size  is 
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60x60.  Because  the  noise  formula  is  valid  only  under 
"steady  state"  conditions,  the  actual  output  is  taken  over 
a  32x32  portion  of  the  output  array,  ignoring  a  band  of 
width  of  14  along  each  of  the  four  edges  of  the  real  output 
array.  The  designed  SVD/SGK  convolution  filter  was 
convolved  with  the  given  input  array  in  fixed-point 
arithmetic.  The  filter  coefficients  were  represented  by 
floating-point  with  36  bits  of  word-length.  The  standard 
deviation  of  the  actual  errors  produced  at  the  output  with 
rounding  to  N  bits  was  measured  and  compared  with  its 
theoretical  estimates  computed  from  the  noise  formula  of 
Eq.  (3-15).  The  system  of  Fig.  5-4  was  used  to  measure  the 
value  of  for  various  word-length  of  storage  0-21.  The 
system  Hro(z^,z2)  was  implemented  with  floating-point 
arithmetic  with  36  bits  of  word-length.  Table  5-1  shows 
the  experimental  results.  There  is  good  agreement  between 
the  predicted  and  measured  values.  This  confirms  the 
validity  of  our  model  and  a  statistical  approach,  to  analyze 
the  roundoff  noise. 

5.2.2  Filter  Coefficient  Quantization  Effect 

In  Chapter  3,  the  quantization  effect  of  the  filter 
coefficients  was  shown  to  be  not  as  severe  as  that  of 
rounding.  Before  we  present  the  experimental  results,  it 
would  be  beneficial  to  discuss  the  error  measurement  of  a 
pair  of  images.  A  true  comparison  between  a  pair  of  images 
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TABLE  5-1 


Standard  deviation  of  outout  noise  caused  by  rounding 
operations  for  a  prototype  filter 


N 

Theory 

Experiment 

8 

0.552x10 

0. 719xl0"2 

10 

-2 

0.138x10 

0. 191xl0-2 

12 

0. 345xl0~3 

0. 465xl9~3 

14 

0. 862xl0~4 

0. 117xl0"4 

16 

0. 216xl0~4 

0. 291xl0"4 

Lowpass  Filter 


N 

Theory 

Experiment 

8 

0. 329xl0-1 

0. 439xl0-1 

10 

0. 82  3xl0~ 2 

0. lllxlO-1 

12 

0. 206xl0-2 

0. 270xl0-2 

14 

0. 514xl0'3 

0. 685xl0-3 

16 

0. 129xl0-3 

0. 173xl0-3 

Bandpass  Filter 
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should  follow  some  objective  criteria.  Tt  is  desirable 
that  the  objective  criteria  be  mathematically  tractible  and 
reasonably  calculable  so  that  they  can  be  used  as  objective 
performance  functions  to  evaluate  an  image  proeossina 
system.  Considerable  attention  has  been  paid  to  rte 
development  of  such  criteria  [5-31.  Unfortunately,  because 
of  the  complexity  of  the  human  visual  system,  there  arc  no 
universally  accepted  criteria  to  measure  image  fidelity. 
Eut,  the  most  commonly  used  quantitative  measure  of  a  pair 
of  image  is  the  normalized  mean  souare  error  Ih’VFFl,  as 
defined  in  Chapter  2  [5-4).  We  shall  use  the  NMFF  as  nur 
objective  criterion  throughout  this  dissertation.  Tatle 
5-2  shews  the  computed  NMSE  between  f  1  oa  t  i  nq -po  i  r.  t 
arithmetic  with  36  bits  of  word-length  and  fixed-poin*- 
arithmetic  with  different  N  and  M  bits  of  word-length.  in 
all  cases,  the  results  obtained  with  P  =16  bits  arc  close 
to  those  with  full  precision.  Tt  is  concluded  that  16  bits 
of  word-length  to  quantize  filter  coefficients  is 
sufficient  without  reducing  filter  performance 
significantly.  In  Chapter  3,  it  was  shown  that  the  storage 
required  for  the  filter  coefficient  is  far  less  than  that 
required  for  the  data.  We  will  then  consider  that  it  is 
more  practical  to  reduce  the  word-length  required  for  the 
data  storage. 
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TABLE  5-2 


Fixed-point  implementation  error  for  various 

word-length 


8 

1 

12 

16 

30 

8 

mm 

1. 4033 

1.3197 

1.3450 

10 

3. 3241 

0.7768 

0.3174 

0. 3082 

12 

3.3706 

0.6223 

0.0843 

0.0741 

14 

3.3700 

0.5995 

0.0221 

0.0194 

16 

3. 3693 

0.5993 

0.0088 

0.0048 

Floating 

3.3690 

0.5990 

0.00765 

0.0000065 

Lowpass  Filter 


M 

N 

8 

j - 

12 

16 

30 

8 

7.9396 

7.3349 

7.1219 

7.1115 

10 

3.9017 

1.8639 

1.8007 

1.7972 

12 

3.5280 

0. 6642 

0.4529 

0.4463 

14 

3.5184 

0.4887 

0.1123 

0.1081 

16 

3.5196 

0.4842 

0.028 

0.028 

Floating 

3.5190 

0. 4841 

i 

0.00435 

0.0000191 

Bandpass  Filter 


5.2.3  Output  Image  Comparison 


In  order  to  evaluate  the  performance  of  the  SVD/SGK 
convolution  more  precisely,  let  us  define  the  following 
NMSE  factors.  Assuming  that  G  and  F  are  output  and  input 
arrays,  respectively,  then  we  shall  use  the  following 
notation 


G  =  F  »»  H 


(5-la) 


— SVD  -  ®®  — SVD 


(5-lb) 


-SVD/SGK  -  *®  — SVD/SGK 


( 5-lc) 


where  H  is  a  prototype  impulse  response,  H  is  the 

~  ~  SVD 

approximation  of  H  by  retaining  a  few  dominant  terms  in  the 

A 

A 

SVD  expansion  of  H,  and  H  is  the  SVD/SGK  convolution 

~  SVD/SGK 

realization  of  H  .In  Eg.  (5-lc),  the  term  on  the  right 
-SVD  * 

has  been  computed  using  fixed-point  arithmetic.  Then  we 
define 


cl~ 


EE 

A _ i 


G(i' j)_GSVD(i'j) 


2  T 


LE  E|G(i'j))| 


1/2 


( 5-2a) 
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Two  errors  are  involved  in  SVD/SGK  convolution  with 
fixed-point  arithmetic:  is  the  error  caused  by  the  FVP 

approximation,  and  c 2  is  the  error  due  to  fixed-point 
arithmetic.  £3  is  the  total  error.  Table  5-3  summarizes 
the  computed  NKSE  with  different  word-length  of  data 
storage.  In  this  experiment,  the  filter  coefficients  were 
quantized  with  16  bits,  and  the  input  was  a  random  array 
with  correlation  coefficients  of  0.95.  Returning  to 
Eq.  (5-2),  we  shall  derive  an  upper  bound  of  the  total 
error  £3.  Since  the  total  error  £3  is  contributed  by 
as  well  as  £2  >  this  bound  will  be  very  useful  in  SVD/SGK 
convolution  implementation.  Let  us  rewrite  Fq.  (5-2)  in 
terms  of  a  matrix  Euclidean  norm,  which  is  defined  to  be 

II  G  ||  =  [j]^|G(i,j)|2  ]V2  (5-3) 

i  j 

Hence,  Eq.  (5-2)  can  be  rewritten  as 


G 


( 5-4  a ) 


/*» 

e3  I!  5  II  2  '  115  -  Ssvd/sgk'I  2  ,5-4c’ 

But  note  that 

If  -"-SVD/SGkH  =  If  -"-SVD+-SVD~-SVD/SGkI'  (5-5) 

Using  the  Schwartz  inequality  on  the  right  hand  side  of 
Eq.  ( 5-5 )  ,  we  have 

If  -“-SVD/SGkH  -  (ll  — SVD  ' I  +ll  -SVD"GSVD/SGkH  )  (5  6) 

Substituting  Eq.  (5-4)  into  Eq.  (5-6)  results  in 

£3  II  G  ||  2<zx\\  G  ||  2+2e1e2  ||  G  1 1 1 1  GgVD||  +  E2  ||  GgVD||  2  <5-7e) 

Therefore 


since 

Returning  to  Table  5-3,  we  can  see  that  the  error  never 
exceeds  the  bound  given  by  Eq.  (5-7b).  However,  the 
fixed-point  implementation  error  e2  could  be  reduced  to 
less  than  1.0  %  NMSE  with  12  bits  word-length  of  storage. 
The  error  is  dominant  in  the  bandpass  filter  case. 
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Obviously,  the  e-^  error  decreases  as  more  terms  are 
retained  in  the  SVD  expansion.  Furthermore,  there  is  no 
reason  to  believe  that  ei  will  be  the  same  order  as  of 
Eg.  (2-19).  For  instance,  of  the  bandpass  filter  is 
0.7825  %  with  a  4-stage  expansion,  Put,  is  3.498  i. 
But,  the  situation  is  quite  opposite  in  the  iowpass  filter 
cases.  Appendix  A  describes  the  relation  between  the 
and  ej.  errors.  As  shown  in  Eq.  (A-16)  of  Appendix  A,  the 
£l  error  is  mainly  attributed  to  the  mean  difference 

A 

between  G  and  G  .  Given  fixed  e,  ,  the  e.  error 

-SVD  k  1 

increases  as  the  mean  difference  increases. 

The  prime  goal  of  this  error  analysis  is  to  reduce  the 

error  and  to  force  the  SVD/SGK  processed  output  closer  to 

the  direct  processed  output.  If  we  correct  the  output  so 

that  m  ~  equals  to  m  .  then  the  e  error  is  the  same 
gSVD  9  1 

order  as  ek.  In  the  following,  we  shall  develop  a  simple 
algorithm  to  force  the  mean  difference  equal  to  zero. 

Assuming  the  filter  is  time- invar iant  and  linear  and 
m^  is  the  mean  of  the  input  array,  then 

mg  =  mf  L2  H(i' j)  (5"8) 

i  j 

But,  the  prototype  impulse  response  matr'x  H  is  normalized 
such  that  LE  H(i,j)  =  1,  therefore, 

i  j 
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Also  , 


=  m, 


*SVD 


2EHsvD(1,j> 


where 


(5-10?) 


K 


— SVD 


-  E  Cmc,  £ 

i=l 


(5-10b) 


Substituting  Ea.  (5-8)  into  Ea.  (5-10)  yields  in  Ea.  (5-11) 


a  =  m. 


U.o-  £  £  HSVDu.j>] 


(5-11) 


where  a  represents  the  mean  difference,  i.e.,  m  -rru 

g  gSVD 

In  order  for  nu  to  be  eaual  to  m„  ,  we  simply  add 
9SVD  '  g 

the  quantity  a  to  every  output  pixel.  This  simple  point 
by  point  operation  will  significantly  reduce  the  error. 

Hence,  the  overall  error  e-j  will  be  reduced. 


Table  5-4  shows  the  effectiveness  of  the  mean 

correction  procedure  in  overall  performance.  Compare  and 

e  ,  where  e  denotes  the  total  error  after  mean  correction. 

4  4  _  _ 

The  extra  computation  of  m  and  2^2-s  H  (i,j)  would  be 

f  i  j  SVD 

fully  justified  in  the  bandpass  filter  case  since  a 
substantial  reduction  in  NKSE  can  be  obtained.  The 
usefulness  of  this  simple  mean  correction  procedure  in  a 
photographic  example  will  be  demonstrated  further  in  the 
next  section. 


104 


TABLE  5-4 


The  NMSE  comparison  of  before  and  after 
mean  correction 


\.  N 

8 

10 

12 

14 

16 

No 

Rounding 

\ 

Before 

1. 3445 

0.3174 

0.0883 

0.0297 

0.0249 

0.0243 

After 

1. 3449 

0.3180 

0.0864 

0.0265 

0.0179 

0.0173 

Lowpass  Filter 


N.  N 

8 

10 

12 

14 

16 

No 

Rounding 

Nv 

Before 

7.9674 

4.0326 

3.5160 

3.4989 

3.4984 

3.4981 

After 

7. 3556 

2.0835 

1.1237 

1.0170 

1.0156 

1.0143 

Bandpass  Filte 


The  experiments  have  been  repeated  with  varying 
correlation  coefficients  of  input  arrays.  The  results  are 
shown  in  Table  5-5.  Experimentally,  it  has  been  concluded 
that  the  fixed-point  implementation  error  is  cuife 
independent  of  the  correlation  coefficient  of  the  input 
array.  These  results  confirm  that  the  model  employed  was 
sufficiently  valid  for  simulation. 

5 . 3  Real  Image  Exper imental  Results 

In  this  section,  photograpic  results,  based  on 
computer  simulation,  for  SVD/SGK  convolution  are  presented. 
The  SVD/SGK  convolution  method  with  a  fixed-point 
arithmetic  has  been  applied  to  the  convolution  of  real 
images  as  a  test  of  its  validity. 

From  previous  experimental  results,  using  random 
number  arrays  as  an  input,  it  was  concluded  that  16  bits  of 
word-length  for  filter  coefficient  quantization  and  12  bits 
for  data  storage,  i.e.,  rounding,  were  sufficient  to  limit 
the  effects  of  quantization  end  roundoff  noise  to  less  than 
1.0  %  NMSE  for  most  practical  cases.  Although  this 
conclusion  is  based  on  the  particular  model  discussed  in 
the  previous  section,  we  shall  use  the  same  word-length  in 


the  experiment 

with 

real 

images . 

Figure  5  -  5a 

shows 

an 

original  aeri 

al  scene 

image.  The 

original  image 

conta 

ins 

256x256  pixels 

wi  th 

each 

pixel  *mpli 

tude  quantized 

over 

the 

integer  range 

0 

to  255.  Tn  th 

e  first  step 

of 

the 
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TABLE  5-5 


Summary  of  experiment  with  vary  inn  correlation 
coefficient  of  the  input  array 


u'N 

0.0 

0.1 

0.  3 

0.  5 

0.7 

0.9 

ei 

0.0807 

0.0771 

0. 0684 

0.0573 

0. 0442 

0.0289 

f  2 

0.1227 

0.1175 

0.1144 

0. 1180 

0. 1072 

e  3 

0.1463 

0. 1391 

0.1369 

0.1304 

0.1273 

0.1187 

e  . 

4 

0.1462 

0. 1393 

0.1273 

J 

0 . 1252 

0.1195 

Lowpass  Filter 


o 

o 

0.1 

0.3 

0.5 

0.7 

0.9 

E1 

2.9375 

2.9332 

2 . 9609 

3.0376 

3.1653 

3.4213 

E2 

0.4248 

0.4201 

0.4288 

0.4369 

0.4485 

0.4483 

£3 

2.9557 

2.9335 

2.9681 

3.050 

3.1625 

3.4440 

G  A 

4 

0.7460 

_ 

0.7562 

0.7561 

— j 

0.7967 

0.8378 

1. 0970 

Bandoass  Filter 
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simulation,  each  pixel  of  the  oriqinal  image  was  normal  j  ?cd 
to  the  range  0.0  to  1.0.  Figures  5-5b  and  5-5c  illustrate 
the  direct  processed  output  with  prototype  lowpess  and 
bandpass  filters,  respectively.  The  direct  processed 
outputs  were  obtained  using  floating-point  arithmetic  with 
36  bits  of  word-length.  A  comparison  of  direct  and  FVP/FCK 
convolution  for  lowpass  and  bandpass  filters  with  N  =  1? 
and  M  =  16  bits,  is  given  in  Figures  5-6  and  c-7, 
respectively.  There  are  no  apparent  differences  in  visual 
results  for  direct  and  FVD/FCK  convolution.  The  measured 
KMSE  and  absolute  difference  image,  multiplied  by  a 
specified  scale  factor,  are  also  presented  to  show  the 
accuracy  of  SVD/SGK  convolution.  In  both  cases,  the 
resulting  errors  are  less  than  1.0  %.  This  experiment 
verifies  the  validity  of  the  model  used  in  the  previous 
section.  Figures  5-8  and  5-9  contain  simulation  results 
for  the  -experiment  of  Figures  5-6  and  5-7  when  the 
word-length  for  data  storage  is  reduced  by  setting  F  =  P . 
Obviously,  the  error  contribution  caused  by  insufficient 
word-length  for  rounding  is  significant. 

To  illustrate  how  the  different  SVD  approximat ions  of 
a  given  prototype  impulse  response  affect  the  outputs. 
Figures  5-10  and  5-11  show  the  SVD/FCK  processed  output 
with  different  K.  In  this  experiment,  K  =  1?  and  F  =  1 6 
were  assumed.  It  is  noted  that  the  filter  with  K  =  1 

corresponds  to  the  FUSE  separable  approximation  of  the 
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Example  of  direct  processing  convolution 

a)  Original 

b)  Lowpass  filter 

c)  bandpass  filter 


Best 

Available 

Copy 


(c) 


(d) 


Figure  5-6.  Comparison  of  direct  and  SVD/SGK 

Convolution  for  lowpass  filter  with 
L=15,  K=3,  M=16  bits  and  N=12  bits. 

a)  Original 

b)  Direct 

c)  SVD/SGK  (NMSE=0. 06398%) 

d)  Absolute  difference  X  scale 
factor  200 
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Figure  5-7.  Comparison  of  direct  and  SVD/SGK 

convolution  for  bandpass  filter  with 
L=ll,  K=4 ,  M=16  bits  and  N=12  bits 

a)  Original 

b)  Direct 

c)  SVD/SGK  (NMSE=0. 8742%) 

d)  Absolute  difference  X  40 
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l 


i 


(c)  (d) 


Figure  5-8.  Comparison  of  direct  and  SVD/SGK 

convolution  for  lowpass  filter  with 
L=15,  K=3,  M=16  bits  and  N=8  bits. 

a)  Original 

b)  Direct 

c)  SVD/SGK  (NMSE=1 . 1037%) 

d)  Absolute  difference  X  200 


Comparison  of  direct  and  SVD/SGK 
convolution  for  bandpass  filter  with 
L=Hf  K=4 ,  M=16  bits  and  N=8  bits. 


a)  Original 

b)  Direct 

c)  SVD/SGK  (NMSE=8. 049%) 

d)  Absolute  difference  X  40 
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Figure  5-10.  Lowpass  SVD/SGK  convolution  with 
K=l,2,  L=15,  M=16  bits  and  N=12 
bits . 


a)  SVD/SGK,  K=1 

b)  Absolute  difference  X  200 

c)  SVD/SGK,  K~2 

d)  Absolute  difference  X  200 
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Bandpass  SVD/SGK  convo 
k-1,2,3,  L=ll,  M=1 6  bi 

a)  SVD/SGK,  K=1 

b)  Absolute  difference 

c)  SVD/SGK,  K=2 

d)  Absolute  difference 

e)  SVD/SGK,  K=3 

f)  Absolute  different 


prototype  impulse  response.  For  the  lowpass  filter,  there 
is  no  significant  visual  difference  among  different  K's. 
Eut,  there  is  a  significant  difference  in  the  bandpass 
filter  . 

Figures  5-12  and  5-13  illustrate  the  effect  of  the 
mean  correction  algorithm.  Although  there  is  an  obvious 
improvement  in  image  quality  in  the  bandpass  filter,  the 
improvement  in  the  lowpass  filter  is  not  noticeable  because 
the  output  mean  before  mean  correction  in  the  lowpass 
filter  is  already  close  to  the  input  mean.  The  measured 
NMSEs  (before  and  after),  computed  means,  and  2Shsvd(1,3) 
are  listed  in  Table  5-6.  1  "l 

For  the  bandpass  filter  with  K  =  1,  before  mean 
correction,  the  SVD  approximation  error  is  so  severe  that 
the  SVD/SGK  processed  output  is  almost  saturated.  After 
mean  correction,  the  output  is  subjectively  satisfying,  and 
the  resulting  NMSF  is  significantly  reduced.  This 
experiment  visually  demonstrates  the  effectiveness  of  the 
mean  correction  procedure. 

5 . 4  Conclusion 

This  chapter  has  presented  experimental  results  of 
SVD/SGK  convolution  using  fixed-point  arithmetic  based  on 
computer  simulation.  First,  the  derived  noise  formula 
predicting  roundoff  noise  has  been  confirmed 
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Comparison  between  before  and  after 
mean  correction  with  L=15,  M=16 
bits  and  N=12  bits  (lowpass) . 

a)  SVD/SGK  (K=l) ,  before 

b)  SVD/SGK  (K=l),  after 

c)  SVD/SGK  (K=2) ,  before 

d)  SVD/SGK  (K=2) ,  after 


Comparison  between  before  and  after 
mean  correction  with  L=ll,  M=16  bits 
and  N=12  bits  (bandpass) . 

a)  SVD/SGK  (K=l) ,  before 

b)  SVD/SGK  (K=l) ,  after 

c)  SVD/SGK  (K=2) ,  before 

d)  SVD/SGK  (K=2 ) ,  after 

e)  SVD/SGK  (K=3) ,  before 

f)  SVD/SGK  (K=3) ,  after 


TABLE  5-6 


Summary  of  experiment  with  real  image 


H 

Before 

After 

£  E  H  ( i  i  ) 

i  j  SVD  ; 

NMSE ( % ) 

Mean 

NMSE (%) 

Mean 

D 

12.6752 

0. 8122 

1.9035 

0.7210 

1.1278 

2 

3.1634 

0.6973 

0. 3658 

0.72006 

0. 9682 

3 

0.0739 

0.7201 

0.0640 

0.72001 

0.9998 

Lowpass  Filter  m^  =  0.7202064 


n 

Before 

After 

i  j  ^SVD^'"^ 

NMSE ( % ) 

Mean 

NMSE (%) 

Mean 

B 

134.6718 

1.7231 

23.3312 

0.7202 

2.3928 

2 

15.5961 

0.6046 

3.3338 

0.7199 

0.8397 

3 

13.8732 

0.6167 

2. 5028 

0.7200 

0.8564 

4 

2.8910 

0.7592 

0.8742 

0.7206 

1.0543 

Bandpass  Filter  m^  =  0.7202064 
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expe  r imenta 

1 1  y  • 

This 

ex  pe  r iment 

verifies  that 

the 

stat ist  ical 

no  i  se 

model 

used  to  analyze 

the  roundoff  error. 

It  has  been 

found 

that  f 

=  16  bits 

and 

N  =  12  bits 

are 

sufficient 

to 

limit 

the  effects 

of 

cuant i ?e t ion 

and 

roundoff  noise  to  less  than  1.0  %  NMSF  in  most  coses.  To 
obtain  a  reasonable  decrease  in  the  NMSF,  a  simple  mean 
correction  algorithm  was  proposed.  The  image  Duality 
impr ovemient  obtainable  by  resetting  rhe  output-  mean  eoual 
to  the  input  mean  has  been  demonstrated.  The  pictorial 
images  resulting  from  FVC/FCK  convolution,  as  shown  in  this 
chapter,  suggest  that  this  technique  may  have  some 


application  in  real-time  image  display  system. 


CKAPTFF  6 


PARAMETRIC  DESIGN  AND  FIXFD-POINT  TNFLFVFNTATI ON 
OF  SVD/SGK  CONVOLUTION  FILTERS 

6 . ]  Int  roduct  ion 

In  this  chapter,  we-  will  consider  the  problem  of 
designing  an  SVD/SGK  convolution  filter  for  which  the 
cutoff  frequency  is  parametrically  variable.  Variable 
cutoff  frequency  filters  have  numerous  applications  in 
image  processing.  For  example,  one  might  sequentially 
obtain  a  best  restored  image  by  changing  the  cutoff 
frequency,  hence  the  frequency  response,  of  the  restoration 
operator . 

Since  filter  coefficients  are  generally  a  function  of 
the  filter  cutoff  frequency,  one  can  change  the  filter 
cutoff  frequency  by  varying  all  of  the  filter  coefficients. 
But  this  procedure  requires  changing  a  number  of 
parameters.  Therefore,  it  is  often  impractical  and  too 
complicated.  It  would  be  more  practical  if  one  could 
construct  a  filter  so  that  the  cutoff  frecuency  is 
controlled  by  only  a  few  parameters,  say  one  or  two. 
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Based  on  the  earlier  work  of  Constand in  ides  [6-1,6-21, 
Schussler  and  Winkelnkemper  [6-3]  were  the  first  to  design 
such  a  variable  cutoff  frequency  digital  filter.  •  jt  has 
been  shown,  that  by  replacing  each  delay  element  in  the 
basic  filter  structure  with  a  first-order  all-pass  network, 
a  transformed  filter  whose  frequency  response  is  identical 
to  that  of  the  basic  filter  on  a  distorted  frequency  scale 
is  obtained.  Unfortunately,  the  method  described  is 
restricted  to  FIR  filters,  and  is  not  applicable  to  TTR 
filters.  Furthermore,  the  resulting  transformed  filter  is 
an  HR  filter  because  by  replacing  the  basic  element,  the 
first-order  all-pass  network  becomes  recursive. 
Consequently,  the  linear  phase  property  of  the  basic  FIR 
filter  is  lost.  But,  the  variation  of  the  cutoff  frequency 
can  be  accomplished.  Cppenheim  et  al .  [6-4]  proposed  a  new 
frequency  transformation  technique  in  which  the  resulting 
transformed  filter  is  still  an  FIR  filter,  and  the  phase  is 
linear  if  the  basic  filter  is  a  linear  phase  FIR  filter. 
By  noting  the  fact  that  the  SVD/SGK  convolution  filter  is 
essentially  a  sum  of  separable  filters,  each  weighted  by  a 
singular  value,  and  each  separable  filter  is  an  outer 
product  of  one-dimensional  column  and  row  convolution 
operators,  it  is  possible  to  extend  the  proposed 
one-dimensional  frequency  tr ans fo rmat ion  technique  to  the 
SVD/SGK  convolution  filters.  We  shall  show  that  this 
approach  is  quite  successful  in  designing  a  variable  cutoff 
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SVD/SGK  convolution  filter. 


In  this  chapter,  we  shall 


discuss  only  a  lowpass-to-lowpass  tr ansformat i on . 

Modification  for  highpass-to-highpass  or 

bandpass-to-bandpass  is  rather  straightforward  in  most 
cases.  We  assume  here  that  the  basic  filter  is  a 
two-dimensional  FIB  filter  with  linear  phase.  The  basic 
concepts  of  frequency  transformation  and  modification  to 
the  SVD/SGK  convolution  filter  are  discussed  in  Section 
6-2.  A  fixed-point  implementation  of  the  variable  cutoff 
convolution  filter  and  experimental  results  are 
described  in  Section  6-3. 

6  •  2  Frequency  Transformation 

of  Linear  Phase  FIR  Filters 

A  one-dimensional  FIB  filter  with  impulse  response  of 
length  2C+1  has  a  frequency  response 


2Q 

h(ei")  =  £  h(m)e'jmu 
m=0 

A  linear  phase  filter  is  symmetrical  so  that 


(6-1 ) 


h(m)  =  h(Q-m)  (6-2) 

for  m  =  0, 1 , .  .  .  ,Q.  Thus 

Q-l 

)  =  e  ^  ®[h(Q)+  Y  2h  (m)  cos  [u  (Q-m)  ]  ] 
m=l 


(6-3) 


Letting  n  =  Q-m,  Eg.  (6-3)  becomes 


Q 

/?(e^U)  -  e  ^]a(n)cosu  (6-4) 

n=0 

where  s(0)  =  h(Q)  and  a(n)  =  2h(C-n)  for  n  =  1,2,...,C. 

Ke  note  that 

Tn (cosu) ■ =  cosnu  (6-5) 

where  T  is  the  n-th  degree  Chebyshev  polynomial  that 
n 

satisfies  the  recursion  formula 

T  .  (x)  =  2xT  (x)  -  T  .  (x)  (6-6) 

n+1  n  n-i 

for  n  =  I,2,...,C.  Thus,  Eq.  (6-4)  can  be  reformulated  as 

Q 

Me^u)  =  e  b  (n)  (cosu)  n  (6-1) 

n=0 

The  new  coefficient  b(n),  for  n  =  0,1,. ..,C>  is  obtained 

from  the  Chebyshev  polynomial  recursion  formula  of 

Eq.  (6-6).  The  basic  approach  [6-5]  to  the  variable  cutoff 

linear  phase  filter  is  to  use  the  tr ansformat ion 

P 

cosu  =  2^Ak(cos8)k  (6-8) 

k=0 

where  u  and  8  are  the  freauency  variables  of  the  basic  and 
transformed  filters,  respectively.  The  transformation 
described  above  preserves  the  frecu^ncy  response  of  the 
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basic  filter,  although  the  freouency  scale  is  distorted  by 
the  transformation.  Ey  substituting  the  transformation  of 
Eg.  (6-8)  into  Eq.  (6-7)  ,  the  frequency  response  of  the 
transformed  filter  is  found  to  be 


(iT(e 


j  (3 , 


e-j6QP 


'  Q  P 

Y.  b(n)  Y  A]:  (cos  6)  k 

Ln=0  k=0 


(6-9) 


From  Eq.  (6-9) ,  it  is  noted  that  the  impulse  response 
dimension  of  the  transformed  filter  is  now  2QP+] .  Ey 
appropriately  controlling  the  parameters  A^,  for 

k  =  0,1,. ..,F,  the  cutoff  frequency  of  the  transformed 
filter  can  be  varied. 

If  P  =  3,  then  Eq.  (6-8)  becomes 

cosu  =  Aq+A^cosB  (6-30) 

and  for  F  =  2,  the  transformation  assumes  the  form 

2 

cosu  =  A0+A^cos8+A2cos  B  (6-13) 

We  shall  call  the  transformations  of  Eq.  (6-30)  and 
Eq.  (6-11)  first-order  and  second-order  transformations, 
respectively.  The  nature  of  the  first-order  transformation 
is  depicted  in  Fig.  6-1.  For  first-order  tr ansformat  ion , 
if  one  is  interested  in  increasing  the  cutoff  frecruency, 
i.e.,  u  <B  ,  where  u_  and  e~  correspond  to  the  cutoff 

C  —  c  U  G 

frequency  of  the  basic  and  transformed  filters, 

respectively,  one  may  prefer  to  constrain  the  transformed 
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>^(u) 


frequency  response  magnitude  at  6  =  0  to  be  equal  to  that 
of  the  basic  filter.  Mathematically,  it  can  be  shown  as 


VejB> 


=  Heju) 


6=0 


u=0 


(6-12) 


in  which  case,  A ^+A ^  =  ].  But,  in  order  to  ensure  that 
|cosu|<l,  Aq  should  lie  in  the  range  of 

0  <  A.  <  1  (6-13) 

—  0 

By  changing  the  frequency  control  parameter  AQ  from  2ero  to 
unity,  we  can  obtain  a  transformed  filter  whose  cutoff 
frequency  is  given  by 


6 

c 


cos 


-1 


cos  yAo 

1“A_ 


(6-14) 


In  other  words,  if  we  wish  to  increase  the  filter  cutoff 

frequency  from  u  to  B  with  u  <  B  »  then  the  control 

c  c  c  —  c 

parameter  AQ  is  obtained  from  Eg.  (6-15)  as 


A 


0 


cosB  -  cos  u 
c _ c 

cosB  -1 
c 


(6-15) 


where  0<  A^<  1.  To  decrease  the  cutoff  frequency  of  the 
basic  filter,  the  core espondence  is 


Equation  (6-16)  leads  to  the  constraint  that  A  ^  =  1-AQ;  the 
parameter  A  is  restricted  to  the  range  of 


■1  <  Aq  <  0 


(6-i7: 


The  resulting  transformed  filter  cutoff  freauency  is  given 
by 


6  =  cos 

c 


-ircos  uc~Ao I 

L  1  +  Aq  J 


(6-18) 


Let  us  associate  the  complex  variable  z  with  the  basic 
filter  system  function  H(z)  and  the  complex  variable  Z  with 
the  transformed  filter  system  function  HT(Z).  Then,  the 
transformation  of  Eq.  (6-8)  is  equivalent  to 

P 


- 1  v^>k 

k=0 


(6-19) 


If  the  filter  is  implemented  as  a  cascade  of  SGK  filters, 
it  is  noted  here  that  the  SGK  filter  should  be  symmetrical 
because  the  transformation  is  applicable  only  to  a  linear 
phase  filter.  In  Chapter  2,  it  was  shown  that  the  complex 
zeros  of  H(z)  should  be  grouped  together  in  conjugate  pairs 
to  ensure  that  all  kernels  for  3x1  SGK  filters  are  real. 
But,  a  resulting  3x1  SGK  filter  may  not  be  linear  phase. 
For  example,  the  3x1  SGK  filter  from  grouping  complex 
conjugate  pair  zeros  not  on  the  unit  circle  will  not  be 
symmetrical.  In  order  to  ensure  that  all  SGK  filters  are 
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linear  phase,  complex  zeros  not  on  the  unit  circle  should 
be  grouped  together  in  groups  of  four,  corresponding  to  the 

complex  conjugates  and  reciprocals,  i.e.  ,  a,  a  ,  I  , .  As 

a  a* 

a  consequence,  H(z)  will  have  fourth-order  SGK  filter  with 
system  function  of  the  form 

r?+i 

H.(z)  =  1-2  (— — )cos0^z  '*'+(r?+  +  4cos20.)z  2 

ri  1  ri 

2  ..  '  (6-20) 
r . +1  . 

_ ,  l  ,  a  -3,-4 
-2 (-p- — )cos0^z  +z 

where  and  0^  are  the  magnitude  and  phase  of  one  of  the 
complex  zeros  not  on  the  unit  circle. 

The  same  rule  of  zero  grouping  in  Chapter  2  can  be 
applied  to  the  real  zeros  and  complex  zeros  on  the  unit 
circle.  Therefore,  we  can  obtain  a  realization  of  H(z)  in 
terms  of  a  cascade  of  second-  or  fourth-order  linear  phase 
SGK  filters.  The  z-trensform  of  the  second-  or 
fourth-order  SGK  filters  can  be  written  as 

2  -1 

Hi(z)  =  X>i(n)  (6-21) 

n=0 

But  H(z)  can  be  factored  in  the  form 


H(z)  =  n  H • ( z ) 
i=l 


Q1  2 

■  i?!  [S*>i 

1  1  n=0 


(n) 


(6-22) 


\ 


i, 


i 
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for  Ct<C.  To  obtain  s  variable  cutoff  linear  phase- 
filter,  based  on  the  SVD/SGK  convolution,  each  FGK  filter 
is  transformed  in  the  manner  described  earlier.  The 
Z-transform  of  the  transformed  filter  is 


HT  (Z) 


n  h  (zi 

i  =  l  1  i 


W1 

n 

i=l 


2 

Eb(n) 

n=0 


P 


E\ 


(6-2?) 


Therefore,  the  coefficients  of  the  transformed  FCK  filter 

are  expressed  in  terms  of  the  parameters  P ^  anc  the 

coefficients  of  the  basic  filter  (see  Appendix  F)  .  Ry 

controlling  the  parameter  P  ,  the  transformed  filter  cutoff 

k 

freauency  can  be  varied.  Before  we  present  the 
experimental  results,  let  us  define  a 


6  -u 

R  =  -~-u  C  X  100 
c 


(6-24) 


which  will  be  used  to  describe  the  degree  of  the 
tr snsf o rma t  ion  .  Figure  6-2  show’s  the  frequency  response  of 
a  typical  lowpass  filter  and  the  parameters  that  define  it. 
The  three  parameters  A-|_  ,  a2  »  and  Af  characterise  the 
frequency  response  of  the  filter.  If  the  parameters  for 
the  transformed  filter  are  close  to  those  of  the  basic 
filter,  then  the  transformation  will  adequately  preserve 
the  frequency  response  of  the  basic  filter. 
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0.2 


stop  bond  rA2 


In  order  to  verify  the  first-order  transformation,  the 
following  experiment  was  performed.  The  basic  filter  is  a 
one-dimensional  linear  phase  lowpass  filter  with  an  impulse 
response  length  of  15.  The  measured  cutoff  frequency  of 
the  basic  filter  is  0.6739.  Throughout  this  chapter,  the 
specified  frequency  is  normalized  to  the  range  of  ( 0 ,  it  ). 
Figure  6-3  shows  the  frequency  responses  of  the  transformed 
filter  with  the  parameter  varied  from  0.1  to  0.8.  Table 
6-1  summarizes  the  filter  parameters  and  the  desired  and 
measured  cutoff  frequencies  of  the  transformed  filters. 
There  is  excellent  agreement  between  the  two  values.  But 
it  is  observed  that  the  first-order  transformation  does  not 
adequately  preserve  the  frequency  response  of  the  basic 
filter  as  AQ  goes  to  1.  In  the  case  of  AQ  >  0.6,  the 
resulting  transformed  filters  can  not  be  considered  to  be 
lowpass  filters,  because  the  first-order  transformation 
does  not  constrain  the  frequency  response  at  u  =  8  =  ti  . 
Experimental  evidence  shows  that  there  is  a  trade  off 
between  R  and  the  preservation  of  the  frequency  response  of 
the  basic  filter.  To  preserve  the  frequency  response  of 
the  basic  filter  more  adequately,  R  should  be  relatively 
small.  Thus,  the  penalty  paid  for  large  R  is  that  the 
transformed  filter  does  not  preserve  the  frequency  response 
of  the  basic  filter  as  shown  in  Fig.  6-3. 

As  alternative  to  the  first-order  transformation,  one 
may  apply  a  second-order  transformat  ion .  Tn  the 
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Figure  6-3.  First-order  transformation 


TABLE  6-1 


List  of  the  transformed  filters  and  their  cutoff 
frequencies  usinc?  first-order  transformation 


Ao 

^2 

Cutoff 

Frequency 

R  ( % ) 

A1 

■  j  L 

Des  i  red 

Measured 

0.0* 

0,0 

0. 1 lxio'2 

0.  367 

0. 6739 

0.6739 

0.0 

0.1 

0.0 

0.  39xl0~ 3 

0.  398 

0.7119 

0.7127 

5.6  1 
| 

0.2 

0.0 

-0 . 22xl0-2 

0.419 

0. 7572 

0.758 3 

12  .  4  j 

0.3 

0.0 

0. 37xl0-2 

0.471 

0.8125 

0.31  37 

20 . 6 

0.4 

0.0 

0.  19x1 0~  3 

0.523 

0.8319 

0.3833 

30  .  9 

0.  5 

0.0 

0. 14xl0_1 

0.701 

0.9730 

0.8571 

44.4 

0.6 

0.0 

0. 38xl0~  1 

0.754 

1. ' 001 

1.  1027 

63.2 

0.7 

0.0 

0.  14x10° 

N.  A. 

1.2960 

1.1299 

92.  3 

0.8 

0.0 

0.38x10° 

N .  A . 

1.6640 

1.6703 

146.9 

0.0  means  a  basic  filter. 


second-order  t  r  ens  f  c  r  m  a  t  ion  ,  there  are  three  parameters  , 


A  ,  A  ,  and  A9  to  be  controlled.  For  the  case  in  which  the 
cutoff  frecuency  of  the  transformed  filter  is  qreater  than 
or  equal  to  the  cutoff  frecuency  of  the  basic  filter,  we 
can  put  another  constraint  on  the  t r an s f o rma t i on  .  That  is 


''T(eJ>  ) 


/  ,  "1  LI  .  i 

=  *i  (e J  )  | 


u=^ 


(6-25) 


Py  imposing  the  constraint  of  Ea.  (6-25),  we  can  preserve 
the  frequency  response  of  the  basic  filter  better  than  when 
the  first-order  transformation  is  used.  Put  we  shall  show 
that  the  second-order  transformation  severely  restricts  the 
range  of  t r ans fo rma t ion  .  By  using  a  similar  analysis,  as 
used  in  the  first-order  transformation,  it  can  be  shown 
that  the  parameter  A Q  is  restricted  to  the  range 


0  s  Ao  s  i 

-i  <  fl0  <  0 


u  <  6 
c  ~  c 


u  >  6 
c  c 


(6-26) 


and  the  desired  cutoff  frequency  BQ  is  given  by 


6  =  cos 

c 


-1 


1-  1-4Aq  (cosuc-Aq) 

2A„ 


(6-27) 


Detailed  derivations  of  Fas.  '6-26)  and  (6-27)  are 

given  in  Appendix  C.  Although  the  t r ans fo rma t i on  is 

achieved  by  varying  the  parameter  A  Q,  the  maximum  or 

minimum  F  can  be  obtained  when  Aq  =  1/2  or  -1/2, 
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respectively.  By  substituting  Aq  =  +1/2  into  Fo.  (6-27), 
we  obtain  the  maximum  or  minimum  attainable  cutoff 
frequency  with  the  second-order  transformation.  The 
relationship  between  uc  and  maximum  or  minimum  attainable 
Q  is  shown  in  Fig.  6-4.  Figure  6-5  shows  the  results  with 
the  second-order  tr ansformat i on .  Table  6-2  shows  the 
measured  filter  parameters.  The  basic  filter  is  the  same 
as  previously  described.  Upon  comparison  of  the 
first-order  and  the  second-order  t r ans fo rma t i ons  ,  the 
second-order  transformat  ion  is  seen  to  adequately  preserve 
the  frequency  response  of  the  basic  filter  if  3^  is  within 
the  transform  range.  Eut  ,  the  resulting  transformed  filter 
has  an  impulse  response  length  of  4Q+1 ,  instead  cf  20+1 
obtained  as  with  the  first-order  tr ansformat ion . 


A  second-order  transformation  is  still  possible  even 
if  the  desired  cutoff  frequency  6C  is  out  of  the  transform 
range.  To  extend  the  transform  range,  for  the  case  of 
uc  <  Bc,  the  constraint  given  in  the  Ea.  (6-25)  is  forced 
to  be  satisfied  at  B  =  ,  where  0  <  <  tt  ,  rather  than 
at  6=  tt.  The  price  paid  for  relaxing  the  constraint  is 
that  the  transformed  filter  characteristics  are  sacrificed 


to  some  degree, 
range  of 
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where  the  matching  point  a ^  is  obtained  from 


=  cos 


-1 


4  (cos  u  -1 ) 
_ c 

(cos6c“l) 2 


+  1 


(6-29) 


But,  it  should  be  noted  here  that  equals  tt  whenever  the 
desired  bc  is  within  the  transform  range.  For  the  case  of 
uc  >^c,  the  constraint  given  in  Eq.  (6-12)  also  can  be 
relaxed  by  locating  the  matching  point  at  p.  =  »  ,  wherp 
0  <«;  <  it.  The  range  of  Ag  can  be  shown  to  be 

cos«2-3 

- ^ - <  AQ  <  0  ( 6- 3 Oa ) 


where 


a 


2 


-1 

cos 


4 (cosuc+l) 

(cos  8  +1 ) 2 
c 


1 


( 6- 3  Ob) 


Figure  6-6  illustrates  the  frequency  response  of  the 
transformed  filter  with  the  (relaxed)  second-order 
transformation.  Table  6-3  lists  the  matching  points,  the 
desired  and  measured  cutoff  frequencies,  and  R.  There  is 
also  good  agreement  between  the  desired  and  measured  cutoff 
frequencies.  As  shown  in  Fig.  6-6,  relaxation  of  the 
constraints  in  the  second-order  transformat  ion  allows  the 
basic  filter  to  transform  in  the  desired  manner,  but  the 
transformation  gradually  degrades  the  freauency  response  of 
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1.0 


normalized  frequency 


Figure  6-5. 


Second-order  transformation  examples 
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TABLE  6-2 


List  of  the  transformed  filters  and  their  cutoff 
frequencies  using  second-order  transformation 


Ao 

A1 

A2 

Af 

Cutoff 

Frequency 

R(%) 

Desired 

Measured 

-0.5 

0.0 

0. llxlO-2 

0.272 

0.4788 

0.4794 

-28.9 

-0.3 

0.0 

0. llxlO-2 

0.304 

0.5362 

0.5371 

-20.4 

-0.  1 

0.0 

0.  llxlO-2 

0.346 

0.6181 

0.6194 

-8.3 

0.0 

0.  0 

0. llxlO-2 

0.367 

0.6739 

0.6739 

0.0 

0.1 

0.0 

0. llxlO-2 

0.398 

0.7444 

0.7463 

10.5 

0.  3 

0.0 

0. llxlO-2 

0.408 

0.9477 

0.9489 

40.6 

0.5 

0.0 

. 

0. llxlO-2 

..  .. 

0.356 

1.2252 

1.2283 

81.8 

NORMALIZED  FREQ. 


Figure  6- 


Second-order  transformation  examples 
with  a  relaxed  constraint 


transformed  filters  and  their  cutoff  frequencies  using 
(relaxed)  second-order  transformation 


the  transformed  filter  as  R  increases.  It  is  believed  that 
with  the  second-order  transformation,  there  should  also  be 
a  trade  off  between  F  and  the  preservation  of  the  freauency 
response  of  the  basic  filter. 

To  extend  the  discussed  frequency  transformation 
technique  to  SVD/SGK  convolution  filters,  it  is  noted  again 
here  that  a  SVD/SGK  convolution  filter  is  a  sum  of 
separable  filters,  and  each  separable  filter  is  decomposed 
into  an  outer  product  of  one-dimensional  convolution 
operators  on  the  columns  and  rows  of  the  input  image.  If 
the  basic  SVD/SGK  convolution  filter  is  a  two-dimensional 
FIR  filter  with  linear  phase,  each  SVD-expanded  separable 
filter  is  also  a  two-dimensional  FIR  filter  wi  h  linear 
phase.  Ke  assume  here  that  the  size  of  the  basic  filter  is 
( 2Q+1 ) x ( 2C+1 ) •  Thus,  each  separable  filter  has  the 
property  that 

Hi(n,m)  =  Hi(Q-n,Q-m),  0  <  m,n<  Q  (6-31) 

Since  H^(n,m)  is  separable,  then 

H^^m)  =  h?(n)hf(m)  (6-32) 

Eut,  from  Eq.  (6-31),  H^(n,m)  is  also  decomposed  into 

H  i  ( n  ,  m )  =  h*T  (Q-n)  h^  (Q-m)  (6-3?) 

Therefore,  the  convolution  operators  on  the  column  and  row, 

h^(n)  and  h^(m) ,  are  also  one-dimensional  linear  phase  FIR 

14  5 


Figure  6-7  illustrates  the  frequency  responses  of  the 
SVD-expanded  separable  filters  on  the  horizontal  axis.  The 
prototype  filter  is  a  two-dimensional  linear  phase  lowpass 
filter,  and  the  SVD/SGK  convolution  filter  was  obtained  by 
truncating  the  SVD  expansion  to  4  terms.  The  first-stage 
separable  filter  corresponding  to  the  largest  singular 
value  shows  almost  the  same  frequency  characteristics  as 
the  basic  filter.  But  the  other  separable  filters 
corresponding  to  the  next  largest  singular  values  no  longer 
are  lowpass  filters. 


But  a  variable  cutoff  SVD/SGK  convolution  filter  is 
obtained  by  simply  transforming  each  of  the  one-dimensional 
convolution  operators  on  the  columns  and  rows  of  the  input 
image.  Specifically,  the  z-transform  of  the  SVD/SGK 
convolution  filter  is  given  by 


K  fQl  Q2 

HSVD/SGK(zl'z2)=^Lf  ^1H«.,i(zl)  ?I1H£,j(z2) 
1=1  L  x  3-1  J 


( 6-3  4a ) 


K  ,Q,  2 


HSVD/SGK(Z1'Z2) 


Sl=l  '  >-n=i  «  o  Z  ,  -1 


( 6-3  4  b) 


z_  +  z?  m" 

.  (m)  (— -2  -  ) 


but  0^,0  <Q.  The  Z-transform  of  the  transformed  SVD/SGK 
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Fa.  (6-34)  by 


filter  is  obtained  by  substituting  f  z  +  z~l)  jn 

2 

the  transformation 


V^1 


.  -1 
z2  +  z2 


Zt  +  Z, 


k=0 
P 


Z'2+ 


(6-^5) 


k=0 

If  the  basic  filter  is  linear  phase  with  cutoff 

frequencies,  and  along  the  horizontal  and  vertical 

C1  L  2 

axes,  respectively,  and  if  one  is  interested  in  changing 


the  cutoff  frequencies  to  g 
transformation  is  performed  a 
one-dimensional  case. 


„  and  g^  ,  basically  the 
C1  c2 


i  t 


v;  3  s 


m 


the 


An  example  of  the  first-order  transformation  on  the 

SVD/SGK  convolution  filter  is  shown  in  Fig.  6-8.  The  basic 

filter  has  quadrilateral  symmetry  with  cutoff  freauencies 

u  =  u  =  0.7097.  A  perspective  view  of  the  frequency 
1  c2 

response  is  shown  in  Fig.  6-S.  Unless  stated  otherwise, 

=  n  =  u  and  e  =8  =6  are  assumed  for  the 

C1  c2  c  ci  c2  c 

tr ansformat ion .  Figure  6-8  shows  the  cross-sec t ional  view 

of  the  frequency  response  on  the  horizontal  axis.  In  the 

case  of  u  <  8  ,  the  transformation  works  quite  adequately, 

C  “*  c 

but  not  for  u  >6  .  Table  6-4  summarizes  the  filter 

c  c 

parameters  and  the  measured  and  desired  cutoff  frecueneies. 

Figure  6-10  shows  another  basic  filter,  which  also  possesses 

quadrilateral  symmetry  with  cutoff  frequency  u  =  1.1266. 
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TABLE  6-4 

List  of  the  transformed  SVD/SGK  convolution  filters  and 
their  cutoff  frequencies  using  first-order  transformation 


Filter 


A1 

A2 

Af 

Cutoff  Frequency 
(Horizontal ) 

R  (%) 

Desired  Measured 

0.25 

0.02 

0.272 

0.7097  0.7097 

0.08 

0.04 

j 

0.188 

0.6500  0.6056 

-8.4 

0.01 

0.08  : 

0.272 

I 

0.5800  0.4684 

|  -18.3 

0.02 

0.04 

0.241 

0.8500  0.8515 

1 

•  19.8 

0.02 

0.03 

0.283 

0.9500  0.9520 

33.9 

*Basic  Filter 


■ 


the  frequency 
ic  filter 
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The  results  of  the  first-  and  second-order  transformations 
are  shown  in  Figures  6-11  and  6-12.  Figure  6-11 
corresponds  to  the  horizontal,  and  Fig.  6-12  corresponds  to 
the  diagonal  direction.  Comparison  of  the  first-order  and 
second-order  transformation  shows  that  the  second-order 
transformation  yields  far  superior  results  (See  Table  6-5). 

Another  interesting  transformation  is  the  case  of 

B  ,  g  .  An  example  of  changing  the  cutoff  freouency  so 
C1  f  2 

that  u  <  8C  u  >  :  ,  is  presented  in  Fig.  6-]?, 

C1  1  c2  c2 

and  the  resulting  filter  perspective  view  of  the  frequency 
response  is  given  in  Fig.  6-14.  In  this  experiment,  the 
second-order  transformation  is  used  with  the  same  basic 
filter  shown  in  Fig.  6-10. 

6  •  ^  Fixed-Point  Implementat ion 

Once  again,  it  is  of  great  interest  to  implement  a 
variable  cutoff  SVD/SGK  convolution  filter  with 
special-purpose  fixed-point  arithmetic  hardware.  Since  the 
transformation  is  mainly  concerned  with  the  cutoff 
frequency  of  the  transformed  filter,  the  effect  of 
fixed-point  implementation  on  the  cutoff  frequency  is 
significant.  Experimental  evidence  shows  that  a  lowpass 
filter  with  filter  coefficients  rounded  to  16  bits  is 
sufficient  for  both  first-  and  second-order  t r ans format i on  . 
The  frequency  responses  with  different  word-length  for 

filter  coefficient  quantization  and  with  no  rounding  are 
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Figure  6-12.  Transformation  on  the  SVD/SGK 
convolution  filter  (diagonal) 


J. 
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Figure  6-13.  Second-order  transformation  for  SQ  ¥ 

(horizontal)  1  * 


plotted  in  Figures  6-15  and  6-16.  In  the  esse  of 

first-order  transformation,  the  response  for  the  8-bit  case 
deviates  from  the  ideal  response  sign i f i cant  1 y  at  the 
beginning  of  the  passband,  whereas  no  visible  errors  are 
seen  anywhere  in  the  stopband.  In  the  case  of  the 
second-order  transformation,  the  response  for  the  12-bit 
case  shows  the-  same  characteristics.  Put  the  responses  for 
16-bit  for  both  first-  and  second-order  transformations  are 
almost  the  same  as  the  ideal.  This  experiment  shows  that, 
for  a  lowpass  filter,  the  cascade  form  is  highly  sensitive 
to  inaccurate  coefficients  in  the  passband,  but  the 
behavior  in  the  stopband  is  much  less  sensitive.  In 
addition,  the  second-order  transformation  requires 
greater  word-length,  to  quantise  the  filter  coefficients 
than  the  first-order  tr ansformat ion  docs.  The  basic  filter 
was  the  same  as  shown  in  Fig.  6-10. 

In  order  to  investigate  the  roundoff  noise  effect  on 

t.he  fixed-point  implementation  of  the  variable  cutoff 

SVD/SGK  convolution  filter,  the  random  number  array  with  a 

size  of  46x46  was  used  again  as  an  input.  The  correlation 

coefficient  was  0.95.  Basically,  the  same  scaling 

procedure  was  used  to  prevent  overflow,  and  the  subopt imal 

ordering  algorithm  of  Chapter  4  to  minimize  the  roundoff 

noise.  Theoretical  estimates  of  the  roundoff  noise 

(standard  deviation),  based  on  the  noise  formula  derived  in 

Chapter  3,  were  computed  and  compared  with  the  measured 
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O'l  m 


second-order  transformation  (horizontal) 


values.  Table  6-6  summarizes  the  results.  Tn  this 
experiment,  we  assumed  that  M  =  16  bits  and  that  only  one 
rounding  operation  is  performed  within  the  SC-K  filter. 
Excellent  agreement  between  the  two  values  is  observed.  Tt 
is  believed,  again,  that  M  =  16  and  N  =  1?  are  sufficient 
to  achieve  less  than  1  %  NMSF  for  both  first-  and 
second -order  transformations.  Surprisingly,  it  is  noted 
that  the  required  word-length  for  the  variable  SVD/SC-K 
convolution  filter  is  the  same  as  required  that  in  the 
SVD/SGK  convolution  filter. 

6 . 4  Conclusion 

In  this  chapter,  attempts  have  been  made  to  develop  a 
design  technique  for  variable  cutoff  SVD/SCK  convolution 
filters.  We  considered  first-  and  second-order 
transformations.  Second-order  transformat i on ,  in  general, 
exhibits  better  results,  but  it  inherently  limits  the  range 
of  transformation.  It  has  been  shown  that  second-order 
transformation  is  still  possible  if  the  specified 
constraint  is  relaxed.  But  the  price  paid  for  relaxing  the 
constraint  is  degradation  of  the  frequency  char acter  ist  ic 
of  the  transformed  filter.  Tn  addition,  the  second-order 
transformat ion  doubles  the  size  of  the  transformed  filter. 
The  problem  of  implementing  the  transformed  filter  in 
fixed-point  arithmetic  was  also  discussed.  For  a  lowpass 
filter,  it  was  shown  that  the  cascade  form  is  highly 
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TABLE  6-6 


Standard  deviation  of  the  transformed  filter  caused 
by  rounaina  operation 


M=  1 6 


N 

Theory 

Experiment 

NMSE (%) 

8 

0. 337xl0-2 

0. 548xl0-2 

2 . 10 

10 

0. 841xl0-3 

0. 1 1 9x  10~  2 

1.21 

12 

0. 211xl0~3 

0.295xl0~3 

0. 47 

14 

0. 526xl0-4 

-  4 

0 . 768x10 

0.  32 

16 

0. 132xl0-4 

0. 203xl0-4 

0.08 

First-Order  Transformation 


M=16 


N 

Theory 

Experiment 

NMSE (%) 

8 

0. 379x10 

0.925x10 

16.23 

10 

0.946x1 0 

0.107x10 

1.14 

12 

0.237x10 

0.314x10 

0. 53 

14 

0.592x10 

0.743x10 

0.2^ 

16 

0. 148x10 

0.305x10 

0.13 

Second-Order  Transformation 
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sensitive  to  inaccurate  coefficients  in  the  passband,  but 
not  in  the  stop  band.  Finally,  it  is  believed  that  M  =  ] 6 
and  N  =  12  are  sufficient  to  obtain  less  than  1  %  K'MfF  in 
most  practical  cases. 
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CHAPTER  7 


SUMMARY  AMD  FUTURE  WORK 

In  this  dissertaion,  attempts  have  been  made  to 
describe  a  novel  architecture  for  performing 
two-dimensional  convolution  with  a  minimum  amount  of 
hardware  using  the  concept  of  sequential  FGK  convolution. 

The  singular  value  decomposition  of  an  impulse 
response  of  a  two-dimensional  FIR  filter  has  proven  to  be 
useful  in  designing  two-dimensional  approximating  FIR 
filters  that  can  be  implemented  as  a  cascade  of  3x1 
convolution  operators.  The  usefulness  of  the  SVD  has  been 
demonstrated  by  noting  a  trade  off  between  approximat ion 
error  and  computational  speed.  The  SVD/SGK  convolution 
approach  is  particularly  attractive  when  one  is  interested 
in  implementing  a  two-dimensional  convolution  with  a 
digital  image  display  system.  An  approach  to 
implementation  for  SVD/SGK  convolution  that  employs  a  small 
set  of  relatively  simple  digital  circuits  has  been 
described.  It  has  been  demonstrated  that  the  statistical 
approach  is  very  useful  to  analyze  effects  of  finite 
word-length  in  representing  filter  coefficients  and  signal 
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magnitudes.  A  theoretical  formula  for  predicting  the  total 
roundoff  noise  has  been  derived  and  confirmed 
exper imentally . 

Two  important  issues  involving  the  implementation  of  a 
digital  filter  as  a  cascade  of  second-order  filters, 
scaling  and  section  ordering  for  SVD/PGK  convolution,  were 
also  considered.  We  have  shown  how  the  algorithm  available 
in  the  domain  of  one-dimensional  signal  processing  can  be 
extended  to  two-dimensional  signal  processing.  One 
interesting  result  is  that  roundoff  noise  can  be  reduced  by 
interlacing  row  and  column  oriented  elementary  second-order 
filters. 

Experimental  results  dealing  with  image  convolution 

show  that  12  bits  are  reouired  for  memo.rv  storage  (the  most 

•  •  • 

expensive  part  in  image  display  systems),  and  16  bits  are 
needed  for  filter  coefficient  quantization  if  one  desires 
to  get  results  indistinguishable  from  the  output  using  full 
precision.  These  features  can  be  reduced  if  one  allows 
some  distortion  in  the  image  outputs.  it  has  been  shown 
that  the  quality  of  an  SVD/SGK  processed  output  is  improved 
by  resetting  the  output  mean  to  be  equal  to  the  input  mean. 
Since  it  is  impractical  to  compute  the  output  mean,  a 
simple  algorithm  for  resetting  the  output  mean  was 
proposed.  The  effectiveness  of  the  proposed  algorithm  has 
been  demonstrated  visually. 
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It  has  been  shown  that  parametric  modification  of  the 
cutoff  frequency  of  a  filter  is  possible  with 
transformation.  Easically,  the  approach  developed  in  the 
one-dimensional  case  by  Gppenhe im  et  al  .  [7-11  was  used.  P 
detailed  analysis  for  first-  and  second-order 
transformations  was  made,  and  several  design  examples  were 
presented.  In  the  first-order  transformation,  the 
transformation  could  not  properly  preserve  the  frequency 
response  of  basic  filter  as  the  degree  of  transformation 
increased.  in  the  second-order  transformation,  due  to 
inherent  characteristics  of  trigonometric  functions  in  the 
transformation,  the  transformation  works  only  in  a  i  mited 
range.  In  other  words,  it  is  impossible  to  change  the 
cutoff  frequency  of  the  basic  filter  arbitrarily.  Put  the 
resulting  transformed  filter  shows  a  frequency  response 
almost  identical  to  the  basic  filter  within  the 
transformation  range.  It  was  found  that,  by  relaxing  the 
specified  constraint,  arbitrary  variation  of  the  cutoff 
frequency  of  the  basic  filter  with  the  second-order 
transformation  is  still  feasible.  But  a  degradation  of  the 
transformed  filter  frequency  response  was  also  observed. 
It  is  believed  that  with  both  transformations,  there  should 
be  a  trade  off  between  the  degree  of  transformation  and  the 
preservation  of  the  frequency  response  of  the  basic  filter. 
Finally,  it  was  found  experimently  that  ]?  bits  for  the 
accumulator  memory  and  16  bits  for  the  filter  coefficients 
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are  also  sufficient  to  limit  Quantization  and  roundoff 
noise  effects  to  less  than  1  %  NMSE  in  both  the  first-  and 
second-order  transformations. 

Several  problems  are  worthy  of  further  investigation. 
If  one  is  particularly  interested  in  implementation  speed, 
the  SGK  convolution  approach  is  always  faster  than  the 
SVD/SGK  convolution  approach.  Comparison  of  the  processing 
frame  cycle  required  for  SGK  and  SVD/SGK  convolution  shows 
that  only  0  frame  cycles  are  needed  with  SGK  convolution, 
while  the  SVD/SGK  convolution  requires  2K0  frame  cycles 
when  the  size  of  impulse  response  is  ( 20+1 ) x ( 2Q+1 )  and  K  is 
the  number  of  singular  values  employed.  Finding  a  simple 
analytical  design  procedure  for  an  SGK  convolution  filter 
is  still  problem.  An  alternative  to  one-dimensional 
SVD/SGK  convolution  is  to  use  two-d imens ional  SVD/SCK 
convolution.  Two-dimensional  SVD/SGK  convolution  reduces 
the  computational  speed  by  a  factor  of  2.  In  this  case, 
the  scaling  procedure  should  be  carefully  chosen  to  use  the 
full  dynamic  range  of  the  given  word-length.  It  is  also 
expected  that  two-dimensional  SVD/SGK  convolution  reauires 
a  greater  word-length. 

In  the  previous  approach  to  the  parametric  design,  we 
used  one-dimensional  methods  to  design  a  variable  cutoff 
SVD/SGK  convolution  filter.  In  addition,  the  basic  filter 
was  restricted  to  be  linear  phase.  A  generalized  design 
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technique  for  the  two-dimensional  variable  cutoff  digital 
filter  would  be  useful.  Finally,  entending  SGK  or  SVD/SGK 
convolution  to  the  recursive  approach  would  also  help  solve 
the  two-dimensional  signal  processing  problem. 
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APPENDIX  A. 


then , 


Suppose  F ( k , £ )  is  an  input  array  and  m  denotes  its 

2  2  ^ 
then  and  can  be  also  expressed  by 

°e  =  EE  EE  lH(i'j)-HSVD(i'i)1A(i-i,'j-j,) 

i 1  j  1  i  j 

* [ H ( i 1 , j ’ ) _HgVD ( i  ’  ,  j  ’  )  ] 

and 

°g  3  LEIE 

i '  j  ’  i  j 

where 

=  EE[F(k'«-mfHF(kti(  £  +  j)-nif] 
k  £ 

Combining  Egs.  ( A— 1 1 )  and  (A-12),  and  substituting 
Eq.  (A-10) ,  yields 


'  2 

- 

C  2 

L-  j' 

2  2 

_i _ i_ 

[H  (i 

» j ) -Hsvd f j ) ^  (i~i 1 f  j,j 

’) 

2  2 

2 

2  H(i, j ) A (i-i • , j-j ’ ) 

i’  j’ 

i 

j 

IH (i ■ , j ■ ) -HSVD(i • , j • )  1 

H (i ’ , j  • ) 
•a 


(A-10) 

mean  , 

( A-l 1 ) 

(A-12) 

(A-l 3 ) 

into 

( A-l 4 ) 
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But  the  first  term  in  Fa.  (A-ll)  is  equivalent 
Thus , 

2  2 

£l  =  ck  *  a 

2  u2  -1 1 

c  M  2 

m 

EE  !G(i'j)'GSVD(i'j)  i2 
i  j _ 

2  M2 

m  M 

1_ 

i  D 

It  is  clear  that  will  decrease  if  one  can  set 
to  zero. 


2 

to  ek  . 


f A-l 5 ) 


( A-l 6 ) 


close 

m 
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APPENDIX  B 


Relation  between  Basic  Filter  Coefficients  and 
Transformed  Filter  Coefficients 


Let  h(n)  for  n  =  0,1,...,2C+1  represents  an  impulse 
response  of  the  basic  filter  and  a(n)  for  n  =  0 , 1 , . . . , ?QP+! 
represents  the  impulse  response  of  the  transformed  filter. 
It  was  shown  that  the  Fourier  transform  of  a  symmetrical 
filter  can  be  expressed  as 

Q 

h{e^u)  =  e  ^uQ  ^b(n)(cosu)n  (E-l) 

n=0 


The  new  coe.f f.iqient  b(nj  is  obtained  from  the  Chebyshev 
polynomial  recursion  formula.  The  relation  between  h(n) 
and  b(n)  is  given  by 


b<or 

1 

o 

H* 

O 

-h(or 

b  (1) 

= 

1  0  1 

h(l) 

b  (2 ) _ 

I 

O 

O 

_ 1 

h  (2) 

( B-2a ) 


for  Q  =  1  and 
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“ 

b  ( 0) 

-1  0  1  0-1 

h(0) 

b  ( 1 ) 

0  10  10 

h  ( 1 ) 

b  (2 ) 

= 

2  0  0  0  2 

h  (2) 

b  (3) 

0  10  10 

h  (3) 

b  (4 ) 

-10  10  -1 

h(4) 

- 

. 

- 

( B  —  2  b ) 


for  Q  =  2.  By  expanding  Eq.  (6-9)  with  F  =  ] ,2,  the 
coefficient  a(n)  can  be  expressed  in  terms  of  b(n)  and  the 
control  parameter  A 


If  P  =  1,  then 


*  a  ( C  r 

1 - 

o 

r'jcM 

O 

1 

■  b  ( 0  )" 

a  ( 1 ) 

= 

-1-  A  l 

2  0  2 

b  (1 ) 

A. 

_a  ( 2 )_ 
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o 

Hi- 

o 

.  b  ( 2  ) . 

for  C  =  1  ana 


a(0) 

-  0 

0 

A1 

0 

0-. 

-  b  (0) -j 

a  ( 1 ) 

0 

A1 

4 

4 

A0A1 

A1 

4 

0 

b  ( 1 ) 
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= 

1 

2 

Ao 

2 

AJ  * 

2 

Ao 

2 

1 

2 

b  (2 ) 

a  (3 ) 

0 

A1 

4 

Vi 

4 

0 

b  ( 3 ) 

-a  (4)- 

-0 

0 

^1 

0 

0- 

-b(4)  -1 

4 


for  C=2,  where  A  +A  =  1  for  u  <6  and  -A  +A  = 

01  c-  c  01 

u  c>  6C- 


( E  -  3  a ) 


(E-3b) 


]  for 
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m6  =  Vl+  T  A1A2 

=  ^0  +  ^2 
m7  2  4 

m8  A0  +  ~2  +  A0A2  +  8  A2 
foe  Q  =  2.  But,  AQ=  —A ^  and  A^=  1  when  P  =  2. 

Therefore,  the  relation  between  a(n)  and  h(n)  can  be 
obtained  directly  by  substituting  Eq.  (B-2)  into  Eq.  (B-3) 
for  first-order  transformation  and  into  Eg.  (B-4)  for 
second-order  transformation. 
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APPENDIX  C 


Derivation  of  Eqs.  (6-26)  and  (6-27) 


The  second-order  transformation  can  be  characterized 


by 


cosu  =  Aq  +  A^cosB  +  A2cos  0 


(C-l) 


By  imposing  the  constraints  at  G  and  tt,  we  obtain 


1  A0  +  A1  +  A2 


-1  -  Aq  -  A1  +  A2 

Equation  (C-2)  immediately  leads  to  the  relations 


(C-2) 


A0  -A2 


A1  =  1 


( C-3a) 
( C-3b) 


Substitution  of  Eq.  (C-3)  into  Eq.  (C-l)  results  in 


cosu  =  Aq  +  cos0  -  AqCos  0 


( C-4 ) 


The  range  of  An ,  satisfying, 


-1  <  f(x)  <  1  (C-5 ) 

2 

will  ensure  that  Icos  ul^l  where  f(x)  =  aq4x~^ox  ' 
x*cos 0  but  -11  x  11.  Note  that  f(x)  is  a  quadratic 
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function  in  x  and  always  passes  through  two  points,  (-1, 
-1 )  and  (1,1).  The  only  case  for  which  the  condition  of 
Eq.  ( C  —  5 )  is  being  satisfied  is  shown  in  Fig.  C-l .  That 
is  equivalent  to  solving 


257  1  1  '  Ao  >  0 


( C-6a) 


it  <  -1  '  \>  <  0 


(C-6b) 


Equation  (C-6)  gives  the  range  of  Aq  such  that 


1  <  A  <  I 

2-0-2 


(C-l ) 


but,  C  _<  P.q<_  i  corresponds  to  the  case  of  uc  <  Bc  and 
<  A.  <0  corresponds  to  the  case  of  u  >  8  , 

2  — *  U  —  C  C 

respectively.  AQ  =  0  means  that  the  transformed  filter  is 
identical  to  the  basic  filter. 

The  relation  between  uc  and  8c  is  obtained  by  solving 
Eq.  (C-7),  which  is  a  quadratic  function  in  cos  8.  That 
is 


cos8c 


1+/1-4Aq (cosuc-A0) 


2A, 


(C-8) 


But  the  plus 
|cos 8  |  <  1  . 


sign  in  Eq.  (C-8)  is 


discarded 


since 


f  ( x) 


Figure  C-l.  Illustration  of  condition 


REFERENCES 


[1-1]  L.B.  Jackson,  J.F.  Kaiser,  and  H.S.  McDonald, 
"An  Approach  to  the  Implementation  of  Digital 
Filers,"  IEEE  Transactions  on  Audio  and 
Electroacoustics,  Vol .  AU-16,  No.  3, 

September  1968,  pp.  413-421. 

[1-2]  W.K.  Pratt,  "An  intelligent  Image  Processing 
Display  Terminal,"  Proc.  SPIE  Tech.  Symp. , 
Vol.  27,  San  Diego,  Calif.,  August  1979. 

[1-3]  H.C.  Andrews,  "Digital  Image  Processing,"  IEFF 
Spectrum,  Vol.  16,  No.  4,  April  1979,  pp.  38-49. 

[1-4]  J.F.  Abramatic  and  O.D.  Faugeras,  "Design  of 
Two-Dimensional  FIR  Filters  from  Small 
Generating  Kernels,"  Proc.  IEEE  Conference  on 
Fattern  Recognition  and  Image  Processing, 
Chicago,  Illinois,  May  1978. 

[1-5]  W.  F.Vii.  Mecklenbrauker  and  R.M.  Mersereau, 

"McClellan  Transformations  for  Two-Dimensional 
Digital  Filtering:  II-lmplementation ,"  IEFE 

Transactions  on  Circuits  and  Systems, 
vol.  CAS-23,  No.  7,  July  1976,  pp.  414-422. 

[1-6]  A . V.  Oppenheim,  W.F.G.  Mecklenbrauker,  and 

R.M.  Mersereau,  "variable  Cutoff  Linear  Phase 
Digital  Filters,"  IEEE  Transactions  on  Circuits 
and  Systems,  Vol.  CAS-23,  No.  4,  April  1976, 
pp.  199-203. 

[2-1]  T.S.  Huang,  Picture  Processing  and  Digital 

Filtering  Topics  '  in  Applied  Physics , 
Spr inger -Verlag  ,  New  Yorlt,  1975. 

[2-2]  M.P.  Ekstrom  and  S.K.  Mitra,  Two-Dimensional 

Digital  Signal  Processing ,  Dowden-Hutchinson  and 
Foss,  Inc.,  1978. 

[2-3]  J.W.  Cooley  and  J.W.  Tukey,  "An  Algorithm  for 

the  Machine  Calculation  of  Fourier  Series," 
Mathematics  of  Computation,  Vol.  19,  April  1965, 
pp.  297-301. 


181 


mm. 


*w 


( 

( 


[2-4]  W.K.  Fratt,  Digital  Image  Processing, 

Wiley-Inter sc ience  ,  New  York,  1978. 

[2-5]  R.M.  Mersereau  and  D.E.  Dudgeon, 

"Two-Dimensional  Digital  Filtering,"  proc.  TEFE, 
Vol.  63,  No. 4,  April  1975,  pp.  610-623. 

[2-6]  B.R.  Hunt,  "Minimizing  the  Computation  Time  by 
Using  the  Techniques  of  Sectioning  for  Digital 
Filtering  of  Pictures,"  IEEE  Transactions  on 
Computer,  Vol.  C-21,  No.  11,  Nov.  1972, 
pp.  1219-1222. 

[2-7]  T.G.  Stockham,  Jr.,  "High  Speed  Convolution  and 
Correlation,"  Proc.  Spring  Joint  Computer 
Conference,  1966,  pp.  229-233. 

[2-8]  J.W.  Cooley,  P.  Lewis,  and  P.D.  Welch,  "The 
Finite  Fourier  Transform,"  IEEE  Transactions  on 
Audio  and  Electroacoustics,  Vol.  AF-17,  No.  2, 
June  1969,  pp.  77-86. 

[2-9]  J.H.  McClellan,  "The  Design  of  Two-Dimensional 
Digital  Filters  by  Transformations,"  7th  Annual 
Princeton  Conference  on  Information  Science  and 
Systems,  March  1973,  pp.  247-251. 

[2-10]  R.M.  Mersereau,  W.F.G.  Mechlenbrauker ,  and 
T.F.  Cuantieri,  Jr.,  "McClellan  Transformations 
for  Two-Dimensional  Digital  Filtering: 

I-Design,"  IEEE  Transactions  on  Circuits  and 
Systems,  Vol .  CAS-23,  No.  7,  July  1976, 
pp.  405-414. 

[2-11]  W.F.W.  Mecklenbr auker  and  R.M.  Mersereau, 

"McClellan  Transformations  for  Two-Dimensional 
Digital  Filtering:  II-lmplementat ion , "  IEEE 


Transactions 

on 

Circuits 

and  Systems, 

Vol . 

CAS-23,  NO 

.  7, 

July  1576,  pp. 

414-422. 

[2-12]  J.H. 

McClellan 

and 

D.S.K.  Chan, 

"A  2-D  FIR 

Filter  Structure  Derived  from  the  Chebyshev 
Recursion,"  IEEE  Transactions  on  Circuits  and 
Systems,  Vol.  CAS-24,  No.  7,  July  1977, 
pp.  372-378. 

[2-13]  J.F.  Abramatic  and  C.D.  Faugeras,  "Design  of 
Two-Dimensional  FIR  Filters  from  "Small" 

Generating  Kernels,"  Proc.  IEEE  Conf.  on 

Pattern  Recognition  and  Image  Processing, 

Chicago,  May  1978. 


182 


[2-14]  J.F.  Abramatic  and  O.C.  Fauger as Subopt imal 
Chebyshev  Design  of  2-D  FIR  Filters  from  "Small" 
Generating  Kernels,"  International  Conf.  on 
Digital  Signal  Processing,  Florence,  Sept.  1978. 

[2-15]  J.F.  Abramatic,  "Image  Filtering  by  Sequential 
Convolution  of  "Small"  Generating  Kernels," 
Proc .  of  the  TEEE  International  Symposium  on 
Circuits  and  Systems,  Tokyo,  Japan,  July  1979. 

[2-16]  G.E.  Goulb  and  C.  Peinsch,  "Singular  Value 

Decomposition  and  Least  Sauare  Solution,"  Numer. 
Math.,  Vol .  14,  1970,  pp.  403-420. 

[2-17]  S.  Treitel  and  J.L.  Shanks,  "The  Design  of 
Multistage  Separable  Planar  Filters,"  IFFF 
Transactions  on  Geosci.  Electron,  Vol.  GE-9, 
No.  1,  Jan.  1971,  pp.  10-27. 

[2-18]  R.E.  Twogood  and  S.K.  Mitra,  "Computer-Aided 
Design  of  Separable  Two-Dimensional  Digital 
Filters,"  IEEE  Transactions  on  Acoust.,  Speech, 
and  Signal  Processing,  Vol.  ASSP-25,  No.  2, 
April  1977,  pp.  165-169. 

[2-19]  B.R.  Hunt,”A  Matrix  Theory  Proof  of  the  Discrete 
Convolution  Theorm,"  IEEE  Transactions  on  Audio 
and  Electroacoustic,  Vol.  AU-19,  No.  4, 
December  1971,  pp.  285-286. 

[2-20]  H.C.  Andrews  and  E.F.  Hunt,  Digital  Image 
Restoration,  Prentice-Hall,  New  Jersey,  1977. 

[2-21]  P.  Lancaster,  Theory  of  Matrices,  Academic 
Press,  New  YorkT  1969. 

[2-22]  N.D.  Mascarenhas  and  Vi.K.  Pratt,  "Digital  Image 
Restoration  Under  a  Regression  Model,"  IEFE 
Transactions  on  Circuits  and  Systems, 
Vol.  CAS-23,  No.  3,  March  1975,  pp.  252-266. 

[2-23]  W.K.  Pratt,  "An  Intelligent  Image  Processing 
Display  Terminal,"  Proc.  SPIE  Tech.  Symp. , 
Vol.  27,  San  Diego,  Calif.,  August  1979. 

[3-1]  L.B.  Jackson,  "Roundoff  Noise  Analysis  for 
Fixed-Point  Digital  Filter  Realized  in  Cascade 
or  Parallel  Form,"  IEEE  Transactions  on  Audio 
and  Electroacoustics,  Vol.  AU-18,  No.  2, 
June  1970,  pp.  107-122. 


183 


[3-2]  L.B.  Jackson,  "On  the  interaction  of  Roundoff 
Noise  and  Dynamic  Ranges  in  Digital  Filters," 
B.S.T.J.,  Vol .  49,  No.  2,  Feb.  1970, 

pp.  159-184. 

[3-3]  B.  Liu,  "Effect  of  Finite  Word  Length  on  the 

Accuracy  of  Digital  Filters  -  A  Review,”  IEEE 
Transactions  on  Circuit  Theory,  Vol.  CT-18, 
No.  6,  Nov.  1971,  pp.  670-677. 

[3-4]  D.S.K.  Chan  and  L.R.  Rabiner,  "Theory  of 

Roundoff  Noise  in  Cascade  Realization  of  finite 
Impulse  Response  Digital  Filters,"  B.S.T.J., 
Vol.  52,  No.  3,  March  1973,  pp.  329-345. 

[3-5]  D.S.K.  Chan  and  L.R.  Rabiner,  "Analysis  of 

Quantization  Errors  in  the  Direct  Form  for 

Finite  Form  of  Finite  Impulse  Response  Digital 
Filters,"  IEEE  Transactions  on  Audio  and 

Electroacoustics,  Vol.  A.U-2]  ,  No.  4, 

August  1973,  pp.  354-366. 

[3-6]  A . V.  Cppenheim  and  R.W.  Schafer,  Digital  Signal 
Processing ,  Prentice-Hall,  New  Jersey,  1975. 

[3-7]  S. A .  Trettler,  Introduction  to  Discrete-Time 

Signal  Processing^  John  Wiley  an3  Sons' ", . New 

York,  1976. 

[3-8]  R.E.  Crochiere,  "Digital  Network  Theory  and  its 
Application  to  the  Analysis  and  Design  of 
Digital  Filter,"  ph.D.  Dissertation,  Dept,  of 
Elec.  Eng.,  M.I.T.,  Cambridge,  Mass.,  1974. 

[3-9]  0.  Hermann  and  H.W.  Schuessler  ,  "On  the  Accuracy 

Problem  in  the  Design  of  New  Recursion  Digital 
Filters,"  Arch.  Elek.  Ubertragung,  Vol.  24, 
Heft  11,  Nov.  1970,  pp.  525-526. 

[4-1]  W.  Schussler,  "On  the  Structure  for  Nonrecursive 
Digital  Filters,"  Arch.  Elek  Ubertragung, 
Vol.  26,  Heft  6,  June  1972,  pp.  255-258. 

[4-2]  D.S.K.  Chan  and  L.R.  Rabiner,  "An  Algorithm  for 
Minimizing  Roundoff  Noise  in  Cascade  Realization 
of  Finite  Impulse  Response  Digital  Filters," 
B.S.T.J.,  vol.  52,  No.  3,  March  1973, 

pp.  347-385. 


184 


[4-3]  L.B.  Jackson,  "Cn  the  Interaction  of  Roundoff 
Noise  and  Dynamic  Range  in  Digital  Filters," 
B.S.T.J.,  Vol .  49,  No.  2,  Feb.  1970, 

pp.  159-184. 

[4-4]  L.R.  Rabiner  and  P.  Gold,  Theory  and  Appl i ca t i on 
of  Digital  S  ignal  Processing ,  Prentice-Fall  ,  New 
Jersey,  1975. 

[5-1]  W.K. Pratt,  G.D.  Faugeras  and  A.  Cagelowics, 
"Visual  Disc r im ina t ion  of  Stochastic  Texture 
Field,"  IEEE  Transactions  on  System,  Nan,  and 
Cybernetics,  Vol.  SMC-8,  No.  11,  November  1979, 
pp.  796-804. 

[5-2]  J.B.  Knowles  and  E.M.  Glcayto,  "Coefficient 
Accuracy  and  Digital  Filter  Response,"  IFEF 
Transactions  on  Circuit  Theory,  Vol.  CT-15, 


No .  1 ,  March 

1968,  pp. 

31-41  . 

[5-3] 

J . F .  Harris, 

Sr  .  , 

"Image  Evaluation  and 

Restoration  ," 

J.  Opt 

.  Soc .  Am . ,  56 ,  No  .  5  , 

May  1566,  pp. 

569-574  . 

[5-4]  h.K.  Pratt,  Digital  Image  Processing, 

ftiley-Inter sc ience  ,  New  York,  1978. 

[6-1]  A . C .  Constantinides  ,  "Frequency  Tr ansformaf ions 
for  Digital  Filters,"  Electron  Letters,  Vol.  ?, 
No.  11,  November  1967,  pp.  487-489. 

[6-2]  A.G.  Constantinides,  "Freouency  Transformations 
for  Digital  Filters,"  Flectron  Letters,  Vol. 

No.  7,  April  1968,  pp.  115-116. 

[6-3]  Vv .  Schuessler  and  W.  Kinkelnkemper  ,  "Variable 
Digital  Filter,"  Arch.  Electr.  Ubertr,  Vol.  24, 
No.  11,  November  1970,  pp.  524-525. 

[6-4]  A  .  V .  Gppenheim,  Vv.F.G.  Keck!  enbrauker  ,  and 

R.M.  Mersereau,  "Variable  Cutoff  Linear  Fhase 
Digital  Filters,"  IEEE  Transactions  on  Circuits 
and  Systems,  Vol.  CAS-23,  No.  4,  April  1976, 
pp.  199-203. 

[7-1]  A .  V.  Cppenheim,  K.F.G.  Necklenbrauker ,  and 

P.M.  Mersereau,  "Variable  Cutoff  Linear  Phase 
Digital  Filters,"  IFEE  Transactions  on  Circuits 
and  Systems,  Vol.  CA.S-22,  No.  4,  April  1  976, 
pp.  1SS-203 . 


185 


